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ABSTRACT 


An  existing  mathematical  model  for  calculating  bottom- 
hole  pressure  in  flowing  single-phase  gas  wells  has  been 
modified.  The  modification  improves  the  various  correla¬ 
tions  employed  in  the  model.  Furthermore,  it  extends  the 
application  to  sour  natural  gas  wells.  The  modified  model 
was  tested  against  field  data.  The  discrepancies  between 
the  modified  model  predictions  and  the  pressure  data  are 
traced  to  a  variety  of  factors. 

Several  side  benefits  have  resulted  from  this  inves¬ 
tigation.  These  are: 

1.  A  new  mathematical  correlation  was  developed  for  the 
ideal  isobaric  heat  capacity,  for  sweet  natural  gases 
over  the  range  32  -  2240°F  and  for  sour  natural  gases 
over  the  range  32  -  300°F. 

2.  A  new  mathematical  generalized  correlation  was  devel¬ 
oped  for  the  heat  capacity  departure  over  a  range  of 
reduced  pressure  of  0  -  15  and  reduced  temperature 

of  1.05  -  3.0. 

3.  Mathematical  expressions  were  developed  for  the  partial 
derivatives  of  reduced  densities  and  compres s i b i 1 i ty 
factors  of  natural  gases  with  respect  to  reduced  pres¬ 
sure  and  temperature. 

4.  A  FORTRAN  subroutine  was  developed  that  is  capable  of 
regenerating  the  Moody  friction  factor  chart  for 
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relative  roughness  up  to  0.05  and  Reynolds 
greater  than  4000. 

A  method  is  proposed  for  the  estimation  of 
roughness  of  tubing.  The  method  is  essentially 
modification  of  the  existing  model. 
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INTRODUCTION 


Flowing  bottom-hole  pressure  is  of  importance  in  study¬ 
ing  reservoir  behavior.  Many  methods  have  been  developed  for 
the  calculation  of  flowing  bottom-hole  pressure  from  surface 
measurements  for  single-phase  gas  flow  (1  ,2,3 ,4,5, 6,7,8) . 
These  methods  use  the  Force-Momentum-Balance  differential 
equation  and  employ  various  simplifying  assumptions. 

In  1964,  Dranchuk  and  Quon  (9)  proposed  a  general  solu¬ 
tion  for  the  equations  describing  steady-state  turbulent 
single-phase  gas  flow  in  circular  conduits.  This  method  has 
an  advantage  over  the  others  in  that  it  permits  the  simultan¬ 
eous  calculation  of  the  pressure  and  temperature  profiles  in 
the  flowing  fluid. 

Inspection  of  the  Dranchuk  and  Quon  method  shows  that 
it  employs  a  constant  step-size  integrating  scheme  and 
approximate  formulae  for  evaluating  the  various  required 
terms  such  as:  the  gas  compressibility  factor,  its  partial 
derivatives  with  respect  to  reduced  pressure  and  temperature, 
the  Prandtl  number,  the  gas-film  heat-transfer  coefficient, 
and  friction  factor.  Furthermore  it  handles  only  the  flow 
of  sweet  natural  gas  and  requires  that  the  gas  composition 
be  known. 

The  objectives  of  this  study  were  to, 

1.  Modify  the  Dranchuk  and  Quon  mathematical  model  so  as  to 
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permit  the  estimation  of  absolute  roughness  when  f 1  owing 
bottom-hole  pressure  is  provided  in  addition  to  various 
surface  measurements. 

2.  Extend  the  applicability  of  both  the  original  and  modi¬ 
fied  models  so  as  to  handle  sour  natural  gas  and  flow 
in  the  region  of  parti al ly  developed  turbulence. 

3.  Study,  test  if  possible,  and  choose  more  sophisticated 
computer  adaptable  methods  for  the  evaluation  of  the  gas 
compressi bi 1 i ty  factor,  its  partial  derivatives,  isobaric 
heat  capacity,  friction  factor,  overall  heat- transfer 
coefficient,  gas-film  heat-transfer  coefficient,  Prandtl 
number,  thermal  conductivity,  and  viscosity. 

4.  Introduce  necessary  modifications  to  both  the  original 
and  modified  models  so  as  to  permit  their  use  to  flow 
problems  where  the  gas  gravity  is  the  only  known  property 
of  the  gas  and  where  gas  composition  is  not  available. 

5.  Compare  the  calculated  bottom-hole  pressure  values  by 
means  of  the  extended  model  with  field  data. 


LITERATURE  REVIEW 


Flowing  Bottom-Hole  Pressure 

A  variety  of  methods  has  been  developed  for  calculat¬ 
ing  flowing  bottom-hole  pressure  in  dry  gas  wells  (10,11,12, 
13,14,15,16,17).  Aziz  (18)  has  reviewed  most  of  these  methods. 
These  methods  utilize  the  Force -Momen tum-Bal ance  differential 
equation  and  employ  simplifying  assumptions  such  as  neglect¬ 
ing  the  kinetic  energy  term,  assuming  constant  fluid  tempera¬ 
ture,  etc.  Dranchuk  and  Quon  (19)  suggested  an  alternative 
approach  by  simultaneously  solving  the  Force -Momen tum-Ba 1 ance 
and  the  Total  Energy  Balance  differential  equations.  Young 
(20)  assumed  a  linear  fluid  temperature  profile  in  order  to 
numerically  integrate  the  Force-Momentum-Balance  differential 
equation.  Recently  Dranchuk  and  McFarland  (21)  have  studied 
Young's  method  and  as  a  result  have  proposed  a  modified  solu¬ 
tion  for  the  same  differential  equation.  The  latter  two 
methods  assume  a  constant  friction  factor  as  suggested  by 
Cullender  and  Smith  (22). 

Compressibility  Factor  of  Natural  Gases 

The  gas  compressibility  factor  is  a  coefficient  which 
corrects  for  the  non-ideal  behaviour  of  real  gases.  This  fac¬ 
tor  is  necessary  in  any  calculation  involving  volumes  of 
gaseous  mixtures.  Standing  and  Katz  (23)  presented  a  graphi- 
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cal  correlation  for  the  compress i bi 1 i ty  factor  of  dry  and 
sweet  natural  gases.  Elfrink,  et  al.  (24)  proposed  another 
graphical  correlation  for  the  same  factor  for  sweet  and  two- 
phase  as  well  as  one-phase  natural  gases.  The  former  corre¬ 
lation  has  been  widely  accepted  by  the  petroleum  industry. 
Accordingly  several  fitting  and  interpolation  techniques 
(25,26,27)  have  sought  a  mathematical  represen  tat i on  for  the 
Standing  and  Katz  correlation. 

In  the  last  35  years  many  equations  of  state  have  been 
developed  to  handle  hydrocarbon  gases  and  their  mixtures 
(28,29,30,31).  Wiehert  and  Aziz  (32,33)  presented  a  brief 
review  of  the  various  equations  of  state  used  for  calculating 
natural  gas  compres s i bi 1 i ty  factors.  Among  these  equations 
those  of  Redl i ch- Kwong  (34)  and  Benedi ct-Webb-Rubi n  (35)  have 
received  wide  interest.  Several  investigators  (36,37,38) 
have  introduced  a  few  modifications  to  the  B-W-R  equation 
to  make  it  applicable  at  very  low  temperatures  and  useful 
in  handling  complex  hydrocarbon  mixtures.  Others  (39,40,41) 
have  fitted  the  generalized  forms  of  various  equations  of 
state  to  the  tabulated  Standing  and  Katz  correlation  (42) 
in  attempts  to  predict  the  volumetric  behavior  of  sweet 
natural  gases. 

In  order  to  predict  the  volumetric  behavior  of  dry 
sour  natural  gases,  Robinson,  et  al .  (43)  developed  a  corr¬ 
ection  factor  correlation  to  be  used  with  the  Standing  and 
Katz  correlation.  Wiehert  and  Aziz  (44,45,46)  presented  an 
alternative  approach  which  employs  modified  pseudocri ti cal s 
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of  sour  natural  gas . 

Isobaric  Heat  Capacity  of  Natural  Gases 

The  isobaric  heat  capacity  of  gases  is  usually  calcu¬ 
lated  by  evaluating  the  ideal  isobaric  heat  capacity  and  the 
heat  capacity  departure  from  the  ideal  state  at  high  pres- 
su  res  . 

Kobe,  et  al .  (47)  and  Thinh,  et  al.  (48)  presented  the 
ideal  isobaric  heat  capacity  of  pure  gases  as  a  function  of 
temperature  by  means  of  third  degree  polynomials.  Kothari 
and  Doraiswamy  (49),  through  dimensional  analysis,  showed 
that  the  ideal  isobaric  heat  capacity  for  pure  gases  could 
be  expressed  as  a  function  of  molecular  weight  and  reduced 
temperature.  On  the  other  hand  the  ideal  isobaric  heat  capa¬ 
city  of  a  gas  mixture  may  be  expressed  by  means  of  the  molal 
average  (50). 

The  heat  capacity  departure  from  the  ideal  state  at  high 

pressures  has  been  studied  by  many  investigators.  Dodge 

C 

(51)  developed  a  generalized  correlation  for  the  ratio  ^  . 

P 

Watson  and  Smith  (52),  Lydersen,  et  al .  (53),  Edmister  (54), 

Weiss  and  Joffe  (55),  and  Sherwood  (56)  presented  either 

graphical  or  tabulated  generalized  correlations  for  the 

difference  (C  -  C°).  Comings  and  Nathan  (57)  showed  that 

C 

the  difference  (C  -  C°)  and  not  the  ratio  is  a  function 

P  P  p 

of  reduced  pressure  and  temperature.  As  a  result,  they 
recommended  the  use  of  the  difference  (C  -  C°)  as  a  genera- 

r  r 

lized  correlation.  Furthermore  they  pointed  out  that  the 
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Watson  and  Smith  generalized  correlation  might  give  errors 
as  high  as  50  per  cent.  The  Edmister  generalized  correlation 
may  be  in  error  by  as  much  as  50  to  100  per  cent  (58). 
Sledjeski  (59),  Seifarth  and  Joffe  (60)  used  the  B-W-R  equa¬ 
tion  of  state  for  predicting  (C  -  C°)  of  propane  and  meth¬ 
ane.  They  reported  that  satisfactory  agreement  with  experi¬ 
mental  data  was  obtained  for  temperatures  greater  than  150°F. 

The  Hydraulic  Friction  Factor 

Under  normal  conditions  of  fluid  flow  in  pipe  (normal 
velocities,  single-phase  flow,  confined  conduits,  no  tempera¬ 
ture  difference  between  the  fluid  and  the  boundary  wall)  the 
friction  factor  is  a  function  of  Reynolds  number  and  rela¬ 
tive  roughness  (61,62).  For  commercial  pipes,  the  friction 
factor  decreases  at  a  decreasing  rate  as  Reynolds  number 

increases,  and  once  the  boundary  Reynolds  number  (R  )  has 

eb 

been  reached,  the  friction  factor  remains  constant  (63,64). 

Several  investi gators  (65,66)  have  developed  different 
mathematical  expressions  for  the  friction  factor  in  both  the 
rough  and  smooth  pipe  regions.  The  friction  factor  in  the 
transition  region  was  controversial  until  Colebrook  (67) 
bridged  the  gap  between  the  smooth  and  rough  pipe  regions 
with  his  particular  transition  law.  Later,  Moody  (68)  con¬ 
structed  a  friction  factor  chart  in  which  he  employed  the 
Von  Karmctn  formulae  for  the  smooth  and  rough  pipe  regions 
and  the  Colebrook  transition  law  for  the  transition  region. 

Cullender  and  Smith  (69)  assumed  that  gas  flow  through 
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tubing  in  gas  wells  is  fully  turbulent.  Consequently  they 
suggested  a  constant  friction  factor  as  a  function  of  tubing 
diameter  only. 

Absolute  Roughness  of  Gas  Well  Tubing 

It  is  well  known  that  absolute  roughness  changes  with 
time  due  to  pipe  corrosion.  Colebrook  (70)  found  that  the 
rate  of  deterioration  of  pipe  depends  on  time,  type  of  flow? 
lining  material,  size  of  pipe,  and  the  character  of  the  flow¬ 
ing  fluid. 

Cullender  and  Binckley  (71)  reported  an  absolute  rough¬ 
ness  of  0.0006  inch  for  gas  well  casing.  Smith,  et  al .  (72), 
from  their  own  experimental  data,  reported  a  value  of  0.00065 
inch  for  the  absolute  roughness  of  gas  well  tubing.  A  range 
of  absolute  roughness  of  0.00055-0.0019  inch  for  experimental 
gas  pipe  lines  has  been  reported  by  others  (73). 

Overall  Heat-Transfer  Coefficient 

Heat- trans fer  between  a  flowing  fluid  in  a  pipe  and  a 
solid  separated  by  one  or  more  solid  layers  can  be  conveni¬ 
ently  handled  by  utilizing  the  overall  heat- trans fer  coeffi¬ 
cient  concept.  Expressions  for  the  overall  heat- trans fer 
coefficient  for  uniform  geometries  may  be  found  in  any  heat- 
transfer  reference  (74,75).  Holst  (76)  studied  this  coeffi¬ 
cient  for  the  case  of  steam  injection  down  a  well  assuming 
simple  and  uniform  geometries.  Dranchuk  and  Quon  (77) 
suggested  an  expression  which  requires  a  tri al -and-error 
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procedure  for  the  same  coefficient  in  a  gas  well  system. 

The  various  expressions  for  the  overall  heat- trans fer 
coefficient  require  knowledge  of  the  gas-film  heat- trans fer 
coefficient.  The  latter  coefficient  has  been  exhaustively 
investigated.  Simple,  semi -emp i ri cal  ,  empirical  and  sophis¬ 
ticated  models  have  been  developed  for  the  estimation  of 
the  gas-film  heat-trans fer  coefficient  (78,79,80,81,82,83). 

Thermal  Conductivity  of  Gas  Mixtures 

Reid  and  Sherwood  (84)  have  reviewed  most  of  the  exist¬ 
ing  methods  for  the  estimation  of  the  thermal  conductivity 
of  pure  gases  and  their  mixtures.  Among  these  methods  that 
of  Misic  and  Thodos  (85,86)  for  pure  gases  may  be  effici¬ 
ently  adapted  to  computer  calculations.  Lindsay  and  Bromley 
(87)  proposed  a  generalized  equation  for  the  thermal 
conductivity  of  gas  mixtures  at  atmospheric  pressure.  Stiel 
and  Thodos  (88)  studied  the  effect  of  pressure  on  the  thermal 
conductivity  of  pure  gases. 

Viscosity  of  Pure  and  Natural  Gases 

The  viscosity  of  pure  gases  and  their  mixtures  has  been 
studied  by  numerous  investigators.  Reid  and  Sherwood  (89) 
gave  an  excellent  review  on  the  subject.  Herning  and  Zipperer 
(90)  suggested  an  empirical  formula  for  estimating  the  vis¬ 
cosity  of  a  gas  mixture  from  those  of  the  individual  compon¬ 
ents.  Bicher  and  Katz  (91)  presented  a  graphical  method  for 
natural  gas  viscosity  as  a  function  of  specific  gravity, 
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temperature,  and  pressure.  Carr,  et  al.  (92)  developed 
another  graphical  correlation,  which  utilized  the  law  of 
corresponding  state.  Lee,  et  al.  (93,94)  studied  and  deve¬ 
loped  empirical  formulae  for  the  viscosity  of  sweet  natural 
gases  as  well  as  pure  hydrocarbons . 


. 


MODEL  DEVELOPMENT 


Four  main  models  were  developed  to  meet  part  of  the 
objectives  of  this  study.  The  objective,  development,  and 
requirements  of  each  main  model  are  as  follows: 


Main  Model  A 

The  objective  of  this  model  is  to  estimate  the  absolute 
roughness  of  tubing. 

The  development  of  this  model  is  presented  in  detail 
in  Appendix  A.  It  can  be  briefly  described  as  follows: 

-  The  Force-Momentum-Balance  differential  equation  for 
a  fluid  element  contained  in  a  pipe  segment  may  be  written  as 


V 


(1) 


where  a  is  the  inverse  of  the  momentum  correction  factor. 

(The  value  of  the  momentum  correction  factor  for 
turbulent  flow  is  closer  to  unity  than  for  laminar 
flow  (95).) 

-  The  Total  Energy  Balance  differential  equation  over 
the  same  pipe  segment  may  be  expressed  as: 


6q 


dW 

J 


ttD2  ,  G  x 

~T~  {iC/ 
w 


dh  + 


<u>d<u> 

ajgc 


(2) 
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-  The  equation  of  state: 


PV  =  ZRT 


(3) 


-  Equations  1,  2,  and  3  were  rearranged  to  yield  the 
following: 


A* 


+  B* 


C* 


(4) 


and 


From  which 


D* 


+  E* 


F* 


(5) 


C*  •  E*  -  B*-F* 
A*  •  E*  -  B**  D* 


(6) 


and 


where 


A* • F*  -  C* • D* 
A*  •  E*  -  B*-D* 


A* 

B* 

C* 

D* 

E* 

and  F* 


=  /cn(Tr»pr>Tc 

=  fon  ^  Tr ’ Pr ,Tc 
■ 

-  fcn( Tr,Pr.Tc 

■  ^VPr’Tc 
=  fon( Tr»Pr.Tc 


Mw 

»Q0 ' 

.D) 

Hw 

’Qo : 

» D ) 

M 

w 

,Q0  > 

» D  ) 

M 

w 

,Q  , 
*  ^0 

» D ) 

Hw 

*  Q  o : 

» D ) 

Hw 

»Q0  ! 

» D  ,U ) 

(7) 


(8) 
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dP  dT 

In  order  to  obtain  a  s  o  1  u  t  i  on,  and  must  be  finite. 
Equations  6,  7 ,  and  8  imply  that 


dP  „ 

<fiT“  =  ^cn(Tr»pr  sTc»pc>Mw»Q0,D*d,'S,U^ 


(9) 


and 


dT, 

dT 


■  ^VPr’Tc’Pc’Mw’VD’'t’S-U> 


(10) 


-  For  a  specific  flow  problem,  that  is,  when  pipe 

dimensions  and  position  in  space  (L*,D,s),  flow  rate  (Q  ), 

and  gas  properties  (T  ,P  ,M  )  are  given;  Equations  9  and  10 

c  c  w 

reduce  to 


dP 

=  fcn{  Tr,Pr  ^,U)  OU 

and 

lr  =  <12> 

-  The  overall  heat- trans fer  coefficient,  U,  as  suggested 
by  Dranchuk  and  Quon  (96)  is: 

1  =  _L 

U  "  ho 


+  R* 


(13) 
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where  hQ  =  gas-film  heat- trans fer  coefficient 
=  fcn(Pr,Re,^) 

=  fon(  Tr.Pr.j)  04) 

Pr  =  Prandtl  number 
=  Reynolds  number 

and  R*  =  average  thermal  resistance  of  the  surround¬ 
ings,  ft . ^  sec.  °R/Btu 

-  Combining  Equations  11,  12,  13,  and  14  yields, 

dP 

^  =  fon( Tr.Pr>f,R*)  (15) 

and 

dT 

Or  =  /cn(Tr,Pr>f-,R*)  06) 

-  It  is  assumed  that  %  and  R*  are  constants  but  have 

d 

unknown  values  which  have  to  be  determined. 

Thus  the  main  model  A,  which  describes  single-phase 
gas  flow  in  pipes,  consists  of  two  non-linear  differential 
equations,  Equations  6  and  7  or  15  and  16,  and  two-parameter 
estimations  namely  ~  and  R*.  This  system  of  differential 
equations  has  to  be  solved  simultaneously  for  the  boundary 
conditions: 


At 

L  =  0  P(0)  =  P, 

and 

T(0) 

At 

L  ■  L*  P(L.j  ■  P2 

and 

T(L*) 
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This  model  requires  algebraic  expressions  for  the  gas 
compressibility  factor,  its  partial  derivatives,  friction 
factor,  the  Prandtl  number,  gas-film  hea t- trans fer  coeffici¬ 
ent,  isobaric  heat  capacity,  viscosity  and  thermal  conducti¬ 
vity. 

Data  required  are:  gas  composition,  flow  rate,  tubing 
head  pressure  and  temperature,  bottom-hole  pressure  and 
temperature,  length  and  inner  diameter  of  tubing,  and  ground 
surface  temperature. 

Main  Model  AA  (Method  2) 

The  objective  of  this  model  is  to  calculate  the  flow¬ 
ing  bottom-hole  pressure. 

This  model  is  the  original  Dranchuk  and  Quon  model. 

It  is  a  special  case  of  main  model  A  where  ^  is  assumed  to 

d 

be  as  suggested  by  Cullender  and  Smith  (97).  Thus, 

d 

this  model  consists  of  two  non-linear  differential  equations. 
Equations  6  and  7,  and  one  parameter  estimation  namely  R*. 
This  system  of  differential  equations  has  to  be  solved 
simultaneously  for  the  boundary  conditions: 

At  L  =  0  T(0)  =  T-j 

At  L  =  L*  ^(L*)  =  ^2  an<^  T(|_*)  ~  ^2 

This  model  requires  algebraic  expressions  for  the  gas 
compressibility  factor,  its  partial  derivatives,  friction 
factor,  the  Prandtl  number,  gas-film  heat- trans fer  coeffici¬ 
ent,  isobaric  heat  capacity,  viscosity  and  thermal  conducti- 


' 
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vity. 

Data  required  are:  gas  composition,  flow  rate,  tubing 
head  pressure  and  temperature,  bottom-hole  temperature, 
length  and  inner  diameter  of  tubing,  and  ground  surface 
temperature . 

Main  Model  B 

The  objective  of  this  model  is  to  estimate  the  abso¬ 
lute  roughness  of  tubing. 

This  model  differs  from  main  model  A  in  that  it  assumes 
the  overall  heat-transfer  coefficient  constant  at  some  aver¬ 
age  value;  hence 


1  =  r**  =  constant  (17) 

Consequently  when  Equations  11,  12,  and  17  are  combined, 

dP 

<nr  =  ^(Vpr>t>R**>  (18> 

and 

dT 

gr1  =  /e«(Tr,Pr,f,R**)  (19) 

The  same  boundary  conditions  as  those  in  main  model  A 
are  employed. 

This  model  requires  algebraic  expressions  for  the  gas 


..  . f-  .  . 


if  mod  .  \  bn&  «S.  <f  >no  My  w  y f ^ .  .da  oC 
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compressibility  factor,  its  partial  derivatives,  friction 
factor,  isobaric  heat  capacity,  and  viscosity. 

Data  required  are:  flow  rate,  tubing  head  pressure 
and  temperature,  bottom-hole  pressure  and  temperature,  length 
and  inner  diameter  of  tubing,  ground  surface  temperature, 
and  either  the  gas  composition  or  the  gas  gravity  along  with 
the  mole  fractions  of  N2,  C02>  and  H^S. 

Main  Model  BB  (Method  2) 

The  objective  of  this  model  is  to  calculate  the  flow¬ 
ing  bottom-hole  pressure. 

This  model  is  a  special  case  of  main  model  B  where 

is  assumed  to  be  .  Thus  this  model  consists  of  two 

d 

non-linear  differential  equations.  Equations  6  and  7,  and 
one  parameter  estimation  namely  R**. 

The  same  boundary  conditions  as  those  in  main  model 
AA  are  employed. 

This  model  requires  algebraic  expressions  for  the  gas 
compressibility  factor,  its  partial  derivatives,  friction 
factor,  isobaric  heat  capacity,  and  viscosity. 

Data  required  are:  flow  rate,  tubing  head  pressure  and 
temperature,  bottom-hole  temperature,  length  and  inner  dia¬ 
meter  of  tubing,  ground  surface  temperature,  and  either  the 
gas  composition  or  the  gas  gravity  along  with  the  mole 
fractions  of  N2,  C02,  and  H2S. 
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Special  Cases  of  Main  Models  AA  and  BB 

In  the  following  special  cases  both  main  models  AA  and 
BB  reduce  to  one  model : 


1 .  Adiabatic  Flow  (Method  1): 

There  is  no  heat- trans fer  between  the  flowing  fluid 
and  the  surrounding;  that  is, 

U  =  0  (20) 


This  method  requires  solving  simultaneously  Equations 
6  and  7  for  the  boundary  conditions. 


At  L  =  0 

and 

At  L  =  L* 


T(0)  T1 
P ( L* )  =  P2 


2 .  Complete  Solution  of  the  Force-Momentum-Balance 

differential  Equation  (Method  3): 

This  method  assumes  a  linear  temperature  profile  in  the 
flowing  fluid.  It  requires  solving  Equation  4  for  the 
boundary  condition: 

L  -  L*  P(Li)  =  P2 

Equation  4  may  be  written  as: 

dP  dT 

^  -  (c*  -  B*  ^)/D* 


(21) 
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where 


f T ( L* ) 

V  |_* 


(22) 


3 .  Method  4: 

This  method  is  a  special  case  of  Method  3  where  the 
friction  factor  is  assumed  constant  as  suggested  by  Cullender 
and  Smith  (98). 

Method  4  is  equivalent  to  that  of  Dranchuk  and 
McFarland  (99). 

The  main  models  require  algebraic  expressions  (models) 
for  the  gas  compres s i bi 1 i ty  factor,  its  partial  derivatives 
with  respect  to  reduced  temperature  and  pressure,  isobaric 
heat  capacity,  friction  factor,  overall  heat- trans fer  co¬ 
efficient,  gas-film  heat-trans fer  coefficient,  the  Prandtl 
number,  thermal  conducti vi ty,  and  viscosity. 

The  various  mathematical  expressions  for  these  terms 
are  the  subject  of  the  rest  of  this  section. 

Compressibility  Factor  of  Natural  Gases 

The  compressibility  factor  of  dry  sweet  natural  gas 
may  be  expressed  analytically  by  many  methods.  In  this  study 
the  reduced  form  of  the  B-W-R  equation  of  state  as  suggested 
by  Purvis  (100)  and  Dranchuk,  et  al.  (101)  was  chosen  to 
represent  the  compressibility  factor,  where 


. 
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ZcPr 


A2  ,  A3 


5  \  /l 


1  pT;  =  1  +  (Ai  +  t;  +  7?)pr +  (A4  +  r;)pr 


+  Vi  p5  +  ^p2{1  +  AgP2r)exP(-A8p2r) 

r 


(23) 


and 


A]  =  0.31506237 
A4  =  0.53530771 
Ay  =  0.68157001 


A2  =  -1.0467099 
A5  =  -0.61232032 
A8  =  0.68446549 


A3  =  -0.57832727 
Ag  =  -0.10488813 


Equation  23  may  be  extended  to  sour  gas  by  applying  modified 
pseudocri ti cal  pressure  and  temperature,  when  calculating 
the  pseudoreduced  pressure  and  temperature,  as  suggested  by 
Wichert  and  Aziz  (102),  where 


e3  *  120 


(YC02  +  V°-9  '  (YC02  +  VH2S^'6 


+  15 


v0 . 5  v4 

yh2s  "  th2s 


(24) 


Tc  '  £3 


and 


P  T  * 
c  c 


Vu  c(l  "  c)c 


h2s 


H2S 7  3 
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Partial  Derivatives  of  Compress ibi 1 1 ty  Factor 

The  partial  derivatives  of  Z  with  respect  to  reduced 
temperature  and  pressure  may  be  calculated  by  employing  a 
numerical  technique.  However  mathematical  expressions  for 
these  derivatives  may  be  derived  from  Equation  23  as  shown 
in  Appendix  B.  These  mathematical  expressions  may  be 
expressed  as: 


(Lz— ) 

l3T  'P 
r  r 


=  Z 


1  1  9pr 

TT  "  p7  (Tr;) 


p 

r  r 


(25) 


and 


(iL_) 

^pA 


=  i 


Pr  pr  Kd?r  Jr 


(26) 


The  partial  derivatives  of  reduced  density  with  respect  to 

3Pr  3P~ 

reduced  temperature  and  pressure,  (jf- )p  and  (td- )t  >  may 

d  r  v r  rr  r 

be  written  as: 


3p 


r 

r  r 


2A 


Pr  +  (ai“t  3^pr  +  A4pr  "  3"  pr^  1+A8pr^ea?p^,A8pr^ 


3  2A7  3 


2 

T„  +  2(AnT.  +  A0  +  +  3(A4Tr  +  A5)Pr  + 


1  r  2  T  2/Kr 
r 


A 


6AkA*p!  +  ~~~o  P*(3  +  3AgP 2  -  2A2p^)^p(-AftP2) 


l5"6Kr  T  2  Kr 
r 


8  r 


(27) 
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and 


9pr 

^9P“^T 
r  r 


Tr  +  2<AlTr  +  A2  +  77)pr  +  3(A4Tr  +  Vpr  + 

r 


6A 


7  2 


2  4 


5A6pr  +  ~T  pr(3  +  3A8pr  '  2A8pr)e*p(-A8pr) 


(28) 


Isobaric  Heat  Capacity  of  Natural  Gases 


The  isobaric  heat  capacity  of  natural  gas  may  be  repre¬ 


sented  by: 


C 


(29) 


where 


isobaric  heat  capacity  of  natural  gases,  Btu/lb 
mole.  °R. 


C°  =  ideal  isobaric  heat  capacity  of  natural  gases, 
^m  Btu/lb  mole.  °R. 


heat  capacity  departure  from  the  ideal  state  of 
natural  gases,  Btu/lb  mole.  °R. 


1 .  Ideal  Isobaric  Heat  Capacity: 

Two  models  may  be  used  to  represent  the  ideal  heat  capa¬ 
city  of  natural  gas. 


a .  Model  C°  1  : 

Pm 

In  this  model  , 
is  , 


is  taken  as  the  molal  average;  that 


recurs  f  i  t  3  -ei  r,<  >< 
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n 


(30) 


where 


(31) 


This  model  requires  that  the  gas  composition  be  known, 
b.  Model  C°  2: 

Pm 

This  model  was  developed  in  this  study  to  represent 
the  ideal  isobaric  heat  capacity  of  sweet  natural  gas  as  a 
function  of  gas  gravity  (g)  and  pseudoreduced  temperature 
(Tr);  that  is, 


A {G)  +  B(G)T 


(32) 


r 


where  the  subscript  m^  c  means  a  mixture  of  hydrocarbons . 
This  equation  can  be  extended  to  sour  natural  gas  by  apply¬ 
ing  corrections  to  account  for  the  presence  of  N2>  C02»  and 
H^S  as  follows: 


C 


o 


FC02  in2  H2S 


(33) 


F 


N 


2 


fcn( Yn  ,T) 


2 


(34) 
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/-(yco2>t) 


(35) 


and 


FH2S  =  ^n(YH2S’T> 
where  T  =  temperature  of  gas,  °R. 


(36) 


2 .  Heat  Capacity  Departure: 

A  heat  capacity  departure  model  was  developed  by  employ¬ 
ing  the  reduced  form  of  the  B-W-R  equation  of  state  and  uti¬ 
lizing  thermodynamic  properties  as  shown  in  Appendix  C.  The 
heat  capacity  departure  may  be  expressed  as: 


C  -C 

pm  % 

R 


-1  - 


6A3  6A,  6A-,  3A 


7  p' 


h  +  +  Pr)exP(‘A8pr) 


T  AT 

'r  m8  r  8  r  r 


2A. 


1  +  (Ar^TT)pr  +  A4pr  -  pr(1+A8pr)exP('A8pP 


-i  2 


r  +  2(AlTr  +  A2  +  7%)pr  +  3(A4Tr  +  A5)pr  + 


6A 


5  ,  A7  2 


5A6pr  +  7T  V3  +  3A8pr 


2A8pr ) ( "A8pr ) 


(37) 


Although  this  equation  was  developed  for  sweet  natural 
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gases,  it  may  be  applied  to  sour  natural  gases  as  will  be 
shown  later. 

The  Hydraulic  Friction  Factor 

The  hydraulic  friction  factor  as  presented  by  Moody 
(103)  in  chart  form  has  been  widely  accepted  and  used  for 
estimating  pressure  drop  due  to  friction  in  both  gas  and  oil 
pipe  lines.  This  chart  was  developed  by  using  the  following 
equations  for  the  various  regions: 

1 .  The  Smooth  Pipe  Region: 


n 


=  2  log  (-§75^-)  ;  £*  <  3 


(38) 


where  f  =  friction  factor. 


Rg  =  Reynolds  number 

Q  G 

orr  W0 


=  20014.855  — -j-  (field  units) 


G  =  gas  gravity,  (air  =  1). 


y  =  gas  viscosity,  c.p. 


Qq  =  gas  flow  rate,  MMSCF/Day. 


£ 


=  absolute  roughness,  inch. 


and 


d 


=  inner  diameter  of  tubing,  inch. 
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2 .  The  Transition  Region: 

This  region  is  sometimes  referred  to  as  the  region  of 
partially  developed  turbulence.  The  friction  factor  in  this 
region  was  represented  by  the  Colebrook  transition  law  (104), 
where 


-  =  -  2 log  (|^4  +  -2— 5-1-)  ;  3  <  e*  <  60  (39) 

/~T~  R  /  f 

e 


3 .  The  Rough  Pipe  Region: 

This  region  is  sometimes  referred  to  as  the  fully  deve¬ 
loped  turbulent  region.  The  Von  Karman  rough  pipe  law  was 
employed  to  represent  the  friction  factor  in  this  region, 
where 


yj-  =  2  log  (!y£)  ;  e*  >  60  (40) 


Rather  than  using  e*  to  express  the  boundary  between 
the  transition  and  rough  pipe  regions  one  may  use  the  bound¬ 
ary  Reynolds  number  (R  )  technique  proposed  by  Moody  (105), 

eb 

where 


3200 

e/d 


(41) 


This  method  is  simpler  and  easier  to  use  than  e*  and  was, 
therefore,  chosen  in  this  study. 

Equation  39  is  too  cumbersome  for  hand  calculation  since 


i  c  i  b9*r  1st  ;  r'  i  -  .-102  2  r 


' 
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it  requires  trial-and-error.  For  this  reason  Cullender  and 
Binckley  (106)  suggest  a  simplified  friction  factor  formula 
as  a  function  of  flow  rate,  viscosity,  gas  gravity,  and  pipe 
diameter.  On  the  other  hand  Cullender  and  Smith  (107)  sug¬ 
gested  two  friction  factor  formulae  for  fully  developed  tur¬ 
bulent  flow  (rough  pipe  region),  both  of  which  are  functions 
of  pipe  diameter  only.  Both  groups  assumed  an  absolute  rough¬ 
ness  of  0.0006  inch  in  the  development  of  their  approximate 
formul ae . 

Moody's  chart  was  employed  (in  this  investigation)  to 
represent  the  hydraulic  friction  factor  for  the  following 
reasons  : 

1.  Main  models  A  and  B  require  that  the  friction  fac¬ 
tor  be  a  function  of  relative  roughness  in  order  to  give  an 
estimate  of  absolute  roughness  of  the  tubing. 

2.  Flow  of  natural  gas  in  gas  wells  is  usually  turbu¬ 
lent  and  falls  in  either  the  transition  or  rough  pipe  region. 
However,  during  "Back  Pressure  Tests",  flow  is  more  likely  to 
be  in  the  transition  region,  especially  at  low  flow  rates. 
Therefore,  the  Cullender  and  Smith  friction  factor  equations 
may  not  hold. 

3.  Moody's  chart  has  been  accepted  as  applicable  to 
pipes  with  relative  roughness  of  0.0  -  0.05  ,  while  the 
other  formulae  are  applicable  only  for  pipes  with  absolute 
roughness  of  0.0006  inch. 


. 
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Overall  Heat-Transfer  Coefficient 

The  heat-transfer  from  a  flowing  fluid  through  pipe  to 
the  surroundings  over  a  pipe  segment  (dL)  may  be  expressed 
as : 


-  6q  =  -  tt  D  dL  U(T$  -  T)  (43) 

For  a  composite  cylindrical  wall,  the  overall  heat- 
transfer  coefficient  (U)  may  be  represented  by  three  models: 


1  .  Model  U-|  : 

This  model  can  be  found  in  any  heat-transfer  reference 
(108).  The  model  that  applies  to  gas  wells  may  be  written  as 


1_ 

U 


D 


n 

r  l 

i  i  =  1 


ln(  D./Dj  ) 


(44) 


2.  Model  U2: 

This  model  was  suggested  by  Dranchuk  and  Quon  (109) 
for  gas  well  systems.  It  assumes  a  constant  thermal  resisti¬ 
vity  of  the  surroundings,  where 


1 

U 


R* 


(45) 


3.  Model  U3: 


This  model  is  proposed  in  this  study  for  gas  well 
systems.  It  assumes  that  the  overall  heat-transfer  coeffici- 


' 
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ent  is  constant  at  some  average  value;  that  is, 

1  =  r**  =  constant  (46) 

Gas-Film  Heat-Transfer  Coefficient 

As  can  be  seen  from  the  previous  section,  models  U-j  and 
require  an  estimate  for  the  gas-film  heat-transfer  co¬ 
efficient,  hQ.  Such  an  estimate,  in  turn,  requires  an  esti¬ 
mate  of  the  Stanton  number. 

Many  investigators  have  developed  various  models  for 
the  Stanton  number,  which  is  related  to  the  gas-film  heat- 
transfer  coefficient  by  (110): 


S 


t 


p<u>c( 


(47) 


where  S^  =  Stanton  number. 

3 

p  =  gas  density,  lbm/ft  . 

<u>  =  mean  gas  velocity,  ft/sec. 
and  Cp  =  isobaric  heat  capacity,  Btu/lbm.°R. 

The  Gowen  and  Smith  (111)  model  is  considered  the  most 
sophisticated  in  addition  to  being  the  best  of  all  in  repro¬ 
ducing  experimental  data.  This  model  may  be  expressed  as: 


X 


^"TTS 


(48) 
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where 


X 


0.5 


155(f) 


0.27 


(49) 


R  =  Reynolds  number. 

c 

PV  =  Prandtl  number, 
r 


f  =  friction  factor. 

k  =  thermal  conductivity,  cal /cm.  sec.  K. 

]i  -  viscosity,  gm/sec.  cm  or  poise. 

Cp  =  isobaric  heat  capacity,  cal/gm  mole.  K. 
Mw  =  molecular  weight,  gm/gm  mole. 


(50) 


Thermal  Conductivity  of  Gas  Mixtures 

Thermal  conductivity  is  required  for  the  calculation 
of  Prandtl  number,  which,  in  turn,  is  necessary  for  the  cal¬ 
culation  of  the  Stanton  number.  Thermal  conductivity  of  gas 
mixtures  may  be  expressed  as: 


k  +  ( k 
m  v  m 


(51) 


where  km  =  thermal  conductivity  of  gas  mixtures. 

k°  =  thermal  conductivity  of  gas  mixtures  at 
m  one  atmosphere. 

(k  -k°)  =  thermal  conductivity  departure  of  gas 
m  m  mixtures. 
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The  Lindsay  and  Bromley  generalized  equation  (112) 
can  be  used  for  the  estimation  of  the  thermal  conductivity 
of  gas  mixtures  at  one  atmosphere,  where 


m 


n 

I 

i  =  l 


k! 

i 


1  + 


Ml 


Y  . 

A  -J- 
Aij  Y. 


(52) 


S 


i 


1.5  T 


b  . 
1 


(54) 


y  s.s, 

i  j 


(55) 


=  the  Sutherland  constant  for  i th  component,  K. 

y!  =  low  pressure  viscosity  of  ith  component,  c.p. 

M  =  molecular  weight  of  ith  component,  gm/gm  mole 
wi 

T,  =  normal  boiling  point  temperature  of  ith  compon- 
D  i  e  n  t ,  K . 


and  k°.  =  thermal  conductivity  of  ith  component  at  one 
1  atmosphere,  cal/sec.  cm.  K. 


The  thermal  conductivity  of  pure  components  at  one 


' 


atmosphere  can  be  represented  as  suggested  by  Misic  and 
Thodos  (113,  114)  where 
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For  Hydrocarbons: 


C* 

P  * 

k?  =  — -r1-  (14.52  T 

10  y .  l 


-  5.14) 


T  >  1 
ri  - 


(56) 


For  Non -hydrocarbons : 
0.75 


C° 

o  P 
k  .  =  — , 

1  1  r.  6 


1 


10  y. 


( 1  9 5 Z  -  31 . 9 4 ) T  +(1  6. 83-82. 5Z  ) 
C  i  ri  ci 


(  1  .  524-2. 8Z  ) 
ci 


1  <  T  <3 
—  r .  — 


and 


1 


2 

3 


(57) 


(58) 


with  T  and  P  in  degrees  Kelvin  and  atmospheres  respect- 
ci  ci 

i ve ly . 

The  thermal  conductivity  departure  can  be  represented 
as  suggested  by  Stiel  and  Thodos  (115),  where 
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14.0  exp ( 0 . 534p  ^ ) 


p  <0.5 
> — 


(59) 


13.1  exp{ 0.67  p  )  -  1  .069 


0.5<pr<2.0  (60) 


c 


2.976  exp{ 1 .155  p  )  -  2.016 


2.0<pr<2.8  (61) 


c 


Viscosity  of  Pure  and  Natural  Gases 

The  viscosity  of  natural  gases  is  required  for  the 
calculation  of  the  Reynolds  and  Prandtl  numbers.  It  can  be 
evaluated  by  means  of  the  Lee,  et  al .  correlation  (116), 
where 


ym  =  K  exp 


i/  =  _ ™ 

*  "  (209  +  19  Mw  +  T) 


(9.4  +  0.02  M  JT1  ,5 

w ' 


(62) 


X  =  3.5  +  +  0.01  M 


w 


Y  =  2.4  -  0.2  X 


where  ym  =  natural  gas  viscosity,  micropoise. 

3 

p  =  natural  gas  density,  gm/cm  . 


DOT  U  3  Z\  r  2r!fip 
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T  =  temperature,  °R. 

and  M  =  average  molecular  weight,  gm/gm  mole, 
w 

They  recommend  that  the  density  be  calculated  using  the 
Standing  and  Katz  compressibility  factor  correlation. 

The  viscosity  of  pure  components  at  low  pressure' is 
needed  for  the  calculation  of  thermal  conductivity  at  low 
pressure.  It  can  be  expressed  by  using  the  modified  rigid 
interacting  sphere  model  (117),  where 


y: 


0.0026693 


7m  r 

/  w. 

T7 

a .  Q, 

1  vi 


(63) 


jji-  =  0.697  (  1  +  0.323  In  T*) 
vi 


(118) 


(64) 


y?  =  viscosity  of  ith  component,  c.p. 


(— )  • 
'  K  i 


T  =  temperature,  K. 

e  =  maximum  energy  of  attraction  for  Lennard- Jones 


o 


poten ti al ,  ergs . 


-16 


K  =  Boltzmann  constant,  1.3805  x  10  ergs/  K. 
and  a  =  collision  diameter  for  Lennard-Jones  potential, 

O 

A. 

e 

Values  of  a  and  (-jr^)  are  available  in  the  literature 
(119,120).  Alternatively,  the  Stiel  and  Thodos  empirical 
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equations  (121)  may  be  used. 

If  T  <  1.5  or  2.0,  Reid  and  Sherwood  (122)  recommend 
ri 

the  modification  of  the  value  of  y?  according  to  the  follow¬ 
ing  relation: 


(T  <2 ) 
r . 


* 

T 


(65) 


where  i s  chosen  arbitrarily  to  correspond  to  that  at 

T  =2. 
r . 


Pseudocri ti cal  Properties  of  Natural  Gases 

The  gas  compressibility  factor,  its  partial  derivatives, 
heat  capacity  and  thermal  conductivity  departure  models  employ 
the  law  of  corresponding  states.  Consequently  an  estimate  of 
the  pseudocri ti cal  properties  is  necessary  to  calculate  the 
reduced  pressure  and  temperature.  These  pseudocri ti cal s  may 
be  taken  as  the  molal  average  as  suggested  by  Kay  (123)  when 
the  gas  composition  is  available.  However  when  the  gas  com¬ 
position  is  not  available  and  where  the  gas  gravity  ( G )  is  the 
only  known  property  of  the  gas,  the  pseudocri ti cal s  may  be  esti¬ 
mated  from  the  Brown,  et  al .  correlation  (124).  Flores  (125) 
presented  the  correlation  in  an  algebraic  form  as 
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T  =  171.137  +  313.725  G  (66) 

p  CH . C . 


P  =  695.100  -  40.000  G  »  G  <  0.85  (67) 

P  CH.C. 


=  704.396  -  51.724  G  ;  G  >  0.85  (68) 


Carr,  et  al  .  (126)  have  introduced  corrections  in  order  to 
make  this  correlation  applicable  to  sour  natural  gases  if 
the  mole  fractions  of  N£»  C02  and  H^S  are  known.  Hence  the 
pseudocri ti cal s  of  sour  natural  gases  may  be  estimated  from 


PTc 


P  c 


-  250  Y 


H.C. 


N, 


80  Y 


CO, 


+  130  Y 


h2s 


(69) 


and 


pPc 


P  c 


H.C. 


-  170  Yn  +  440  Y c0  +  600  Y 


H2S 


(70) 


o  r  i  o  b 


PRELIMINARY  MODEL  TESTING 


The  various  algebraic  expressions  selected  in  the  pre¬ 
vious  section  to  represent  compress ib i 1 i ty  factor,  its  par¬ 
tial  derivatives,  isobaric  heat  capacity,  friction  factor, 
and  viscosity  were  tested  by  comparing  with  experimental 
data  or  correlations  available  in  the  literature.  The 
results  of  these  tests  are  presented  in  this  section. 

Compressibility  Factor  of  Natural  Gases 

Table  1  presents  a  comparison  of  the  compressibility 
factor  model  as  expressed  by  Equation  23  with  1350  experi¬ 
mental  data  points  (127,128,129).  These  experimental  data 
are  for  sweet  as  well  as  sour  natural  gases  with  up  to  73.85% 
H2S  and  54.46%  CO^ 

Partial  Derivatives  of  Compressibility  Factor 

The  partial  derivatives  of  the  compressibility  factor 
with  respect  to  reduced  pressure  and  temperature  were  com¬ 
pared  to  values  determined  graphically.  This  was  achieved 
by  using  experimental  compressibi 1 i ty  factors  (130,131). 
First,  a  few  isobars  and  isotherms  were  plotted,  then  the 
slopes  of  these  curves  at  specific  Pr  and  T^  were  calculated. 
Finally  these  results  were  compared  to  the  partial  deriva¬ 
tives  estimated  by  the  algebraic  expressions.  It  was  found 
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that  Equation  25  gives  an  average  absolute  deviation  of  5.84 
per  cent  when  compared  with  5  points,  while  Equation  26 
deviates  by  1.96  per  cent  when  compared  with  16  points. 

Isobar ic  Heat  Capacity  of  Natural  Gases 

1  .  Ideal  Isobar ic  Heat  Capacity 

The  isobar ic  heat  capacity  of  sweet  natural  gases 
which  have  no  nitrogen  may  be  represented  by  Equation  32  as 
suggested  in  Model  Development.  This  is  an  equation  of  a 
family  of  straight  lines,  the  slopes  and  intercepts  of  which 
are  functions  of  the  gas  gravity. 

To  test  the  model  represented  by  Equation  32,  the  com 
positions  of  72  gases  were  collected  from  the  literature  and 
the  equivalent  pure  hydrocarbon  compositions  were  calculated 
Equations  30  and  31  were  then  employed  to  calculate  the 
ideal  isobaric  heat  capacities  over  a  temperature  range  of 
30  to  300°F  at  10°F  intervals,  thus  producing  2016  data 
points.  Data  points  for  each  gas  mixture  were  plotted  as 
shown  in  Figure  1.  Each  gas  mixture  showed  a  straight  line 
relationship  expressed  by  Equation  32: 

Cpm  =  A(G)  +  B(C)T  (32) 

mH.C.  r 

A {g)  and  B(&)  were  found  to  be  adequately  represented  by 
means  of  second  degree  polynomials  in  G  as  shown  in  Fig¬ 
ure  2;  that  is. 
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FIGURE  1.  IDEAL  ISOBARIC  HEAT  CAPACITY 
OF  NATURAL  GASES 
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GAS  GRAVITY ,  G 


FIGURE  2.  THE  SLOPES  AND  INTERCEPTS  IN  EQUATION  32 
AS  FUNCTIONS  OF  GAS  GRAVITY. 


A (G),  (BTU/lb  mole.  °R) 
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A(gO  =  B-]  +  B  2G  +  ^3^^ 


(71) 


and 


B(g)  -  B4  +  B5g  +  BgS2 


(72) 


The  coefficients  in  Equations  71  and  72  were  obtained  by 
minimizing  the  sum  of  squares  of  errors  in  A {g)  and  B (g) 
respecti vely .  These  coefficients  (will  be  referred  to  as 
the  set  of  coefficients  No.  0)  are: 

B]  =  5.596695  B2  =  -2.233480  B3  =  0.807265 
B4  =  -1.003900  B5  =  3.141600  Bg  =  5.758700 

Equations  32,  71  and  72  give  an  average  absolute  devia¬ 
tion  of  0.079  per  cent  and  a  standard  deviation  of  0.00119 
with  respect  to  the  2016  data  points  used  to  determine  the 
coef fi ci en ts . 

Studying  the  functional  forms  of  Equations  34  through 
36  individually  showed  that  they  can  be  expressed  as: 


F 


CO 


1  +  VC02  (By  +  BgT) 


(73) 


2 


F 


1  +  YH2S  (B9  +  B10T> 


(74) 


and 


F 


N 


2 


where  T  «  temperature,  °R. 
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The  coefficients  through  B-j^  were  obtained  by  mini¬ 
mizing  the  sum  of  squares  of  errors  in  C°  (as  expressed 

^m 

by  Equation  33)  and  were  found  to  be: 

B7  =  0.4258600  Bg  =  1.24323  X  10"3 

B]0  =  1.00889  x  10"3  Bn  =  0.3623622 

B13  =  0.97570  X  10‘3  B14  =  2.708217  X  10"3 

Equations  33,  73,  74,  and  75  give  an  average  absolute 
deviation  of  0.629  per  cent  and  standard  deviation  of  0.00989 
with  respect  to  2016  data  points  used  to  determine  the  co¬ 
efficients  B 7  through  B^. 

Equations  32,  71,  and  72  are  applicable  over  a  tempera¬ 
ture  range  of  30-300°F  to  sweet  natural  gas  which  has  no 
nitrogen.  The  same  equations  can  be  used  over  a  broader 
temperature  range;  however,  it  is  recommended  that  other  sets 
of  coefficients  depending  on  temperature  range  of  operation 
be  used.  The  quality  of  fit  of  the  new  sets  of  coefficients 
is  shown  in  Table  2. 

2 .  Heat  Capacity  Departure: 

The  heat  capacity  departure  model  as  expressed  by  Equa¬ 
tion  37  was  compared  with  the  tabulated  Weiss  and  Joffe 
generalized  correlation  (132).  The  model  gives  an  average 
absolute  deviation  of  5.91  per  cent  when  compared  with  150 
points. 


Bg  =  -0.0405600 
B12  =  -0.4661581 
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Tables  3  and  4  present  comparisons  of  the  various  heat 
capacity  departure  correlations  and  the  proposed  generalized 
correlation  as  expressed  by  Equation  37  with  the  experiment¬ 
al  data  of  Workman  (133),  Krase  and  Mackey  (134).  Table  5 
presents  comparisons  of  Equation  37  and  the  various  (C^  -  C°) 
generalized  correlations  with  that  of  Weiss  and  Joffe. 

Table  6  present  comparisons  of  Equation  37  with  the  experi¬ 
mental  data  of  Sage,  et  al  .  (135)  for  two  natural  gases.  It 
also  presents  the  effect  of  any  errors  in  (Cp  -  C°)  on  the 

isobaric  heat  capacity  (Cp). 

The  Hydraulic  Friction  Factor 

Table  7  presents  comparisons  between  the  values  of 

friction  factors  as  suggested  by  Moody  (136)  and  Cullender 
and  Smith  (137).  This  comparison  is  carried  out  for  the 
same  12  test-points  used  to  test  main  models  AA  and  BB. 

Viscosity  of  Natural  Gases 

The  Lee,  et  al .  viscosity  correlation  (138)  was  tested 
with  their  111  viscosity  data  points.  This  correlation  gives 
an  average  absolute  deviation  of  3.71  per  cent  from  experi¬ 
mental  data  and  a  standard  deviation  of  0.052527.  The  maxi¬ 
mum  deviation  encountered  is  16.38  per  cent. 
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TABLE  5 


COMPARISON  OF  EQUATION  37  AND  OTHER  (C  - 
CORRELATIONS  WITH  THE  WEISS  AND  JOFFE 


C°)  GENERALIZED 
CORRELATION 


%  Deviation  from  the  Weiss  and  Joffe  Values 


Correl ati on 

Tr 

=  1.0 

1  .2 

1  .« 

2.5 

Hougen  and  Watson 

59 

to 

62 

26 

to 

55 

33 

to 

64 

43 

to 

63 

Edmi s ter 

1  1 

to 

37 

-1 

to 

-11 

-4 

to 

-8 

-22 

to 

-23 

Lydersen ,  et  a l  . 

-9 

to 

35 

5 

to 

29 

-61 

to 

20 

-30 

to 

-108 

Sherwood 

13 

to 

35 

1 

to 

35 

-25 

to 

5 

8 

to 

66 

Equation  37 

* 

-6 

to 

6 

-17 

to 

0 

-4 

to 

-13 

*  Outside  the  range  of  correlation. 
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where  W  Indicates  well  head  conditions 

S  indicates  bottom-hole  conditions 
+  indicates  absolute  roughness  of  0.0006  inch. 


. 

DATA  AND  RESULTS 


Back  Pressure  Data 

Back  pressure  data  were  used  to  test  main  models  AA 
and  BB.  The  test  consists  of  12  test-points  taken  from  four 
back  pressure  tests  for  three  natural  gas  wells.  Table  8 
presents  a  summary  of  these  data.  The  pressure  and  tempera¬ 
ture  data  are  shown  in  Table  9.  Appendix  F  is  a  tabulation 
of  the  source  data  and  gas  properties  upon  which  this  study 
was  based. 

Accuracy  of  Field  Data 

1  .  Flow  Rate: 

The  flow  rate  was  measured  using  either  critical  flow 
provers  or  orifice  meters,  and  was  reported  at  the  standard 
conditions  of  14.65  psia  and  60°F.  The  accuracy  of  the  flow 
measuring  devices  is  not  available. 

2 .  Length  of  Tubing: 

The  reported  values  of  length  of  tubing  are  those  from 
Kelly  Bushing  (K.B.)  to  the  mid-point  of  perforations 
(M.P.P.).  The  depth  at  which  the  pressure  measuring  device 
was  located  is  not  stated.  The  length  of  tubing  of  interest 
in  this  study  is  that  between  the  well  head  and  the  bottom- 
hole  pressure  measuring  device.  This  may  be  different  from 


50 


51 


TABLE  8 

SUMMARY  OF  FIELD  DATA 


Length  of  Flow  Strings  ft. 

4753  - 

11029 

Flow  Rate,  MMSCF/Day 

0.833  - 

17.359 

Gas/Condensate  Ratio,  MMSCF/bbl. 

0.208  - 

00 

Tubing-Head  Temperature,  °F 

47 

170 

Bottom-Hole  Temperature,  °F 

115 

246 

Tubing-Head  Pressure,  PSIA 

985 

3250 

Bottom-Hole  Pressure,  PSIA 

1104  - 

4250 

Water/Gas  Ratio,  bbl/MMSCF 

0.0 

5.2 

Yco  ,  %  mole 

0.56  - 

2.88 

Yh  s,  %  mole 

0.0 

18.89 

Ym  ,  %  mole 

1.05  - 

4.56 

PRESSURE  AND  TEMPERATURE  DATA 
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the  values  reported  in  Appendix  F  for  length  of  tubing. 

3 .  Pressure  Data: 

The  tubing  head  pressure  was  measured  with  a  dead 
weight  gauge  or  a  Barton  gauge.  The  bottom-hole  pressure 
was  measured  with  an  Amerada,  tester  type  gauge,  or  the 
Sperry  Sun  gauge.  The  accuracy  of  these  pressure  devices 
is  shown  in  Table  10  (139).  The  accuracy  of  the  pressure 
data,  is  presented  in  Table  11. 


TABLE  10 

ACCURACY  OF  PRESSURE  MEASURING  DEVICES 


Measuring  Devices 


Accuracy,  %  of  Rated  Pressure 


Dead  Weight  Gauge 

Barton  Gauge 

Sperry  Sun  Gauge 

Amerada  or  Tester  Type  Gauge 


+  0.0 

+0.5*  to  £  0. 75** 
+0.1 

+0.2*  to  +  0 . 25** 


*  For  rated  pressure  less  than  or  equal  to  2000  psi. 
**  For  rated  pressure  greater  than  2000  psi. 


Results 

Computer  FORTRAN  programs  were  written  for  main  models 
AA  and  BB  (method  2)  and  methods  1,  3,  and  4. 

Table  12  presents  comparisons  between  the  measured  and 


ACCURACY  OF  PRESSURE  DATA 
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calculated  flowing  bottom-hole  pressures  for  all  test-points. 
The  test-points  are  presented  in  order  of  increasing  flow 
rate.  Table  13  presents  the  deviations  of  the  calculated 
from  the  measured  flowing  bottom-hole  pressures,  while  Table 
14  exhibits  the  per  cent  deviations  of  the  calculated  from 
the  measured  pressure  drops.  Table  15  presents  the  various 
flowing  bottom-hole  pressures  calculated  by  main  models  AA 
and  BB. 


COMPARISONS  BETWEEN  THE  MEASURED  AND  THE  CALCULATED  FLOWING  BOTTOM-HOLE  PRESSURES 
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COMPARISON  OF  CALCULATED  BOTTOM-HOLE  PRESSURES  AND  COMPUTATIONAL  COSTS 

FOR  MAIN  MODELS  AA  AND  BB 
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The  gas  composition  was  supplied  to  the  program  as  data. 


DISCUSSION 


Compressibility  Factor  of  Natural  Gases 

The  natural  gas  compressibility  factor  model  as 
expressed  by  Equation  23,  which  is  an  equation  of  state,  can 
predict  Z  with  an  average  deviation  of  1.013  per  cent  from 
experimental  data  as  shown  in  Table  1.  Such  accuracy  is 
considered  acceptable  for  engineering  applications. 

Partial  Derivatives  of  Compressibility  Factor 

The  partial  derivatives  of  the  gas  compressibility 
factor  may  be  calculated  by  employing  a  numerical  scheme  or 
by  using  Equations  25  and  26.  However,  a  numerical  technique 
is  cumbersome  and  costly.  These  equations  are  expected  to 
work  well  since  they  were  derived  from  an  equation  of  state. 
Equation  23.  Robinson  (140)  states  that  the  B-W-R  equation 
of  state  gives  good  predictions  for  the  partial  derivatives 
of  compressibility  factor. 

Isobaric  Heat  Capacity  of  Natural  Gases 

1 .  Ideal  Isobaric  Heat  Capacity: 

The  proposed  model  for  ideal  isobaric  heat  capacity 
of  natural  gases  as  expressed  by  Equations  32,  71,  and  72 
gives  excellent  results  as  was  shown  earlier  in  Preliminary 
Model  Testing.  However,  it  is  limited  in  its  application  to 
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a  specific  temperature  range  and  natural  gases  which  consist 
of  pure  hydrocarbons .  Furthermore,  it  requires  that  the 
hydrocarbon  gas  mixtures  be  rich  in  methane. 

The  ideal  isobaric  heat  capacity  of  sweet  natural  gas 
may  be  expressed  by  the  same  equations  over  a  wider  range 
of  temperature;  however,  new  sets  of  coefficients,  B-j 
through  Bg,  should  be  used.  The  set  of  coefficients  No.  1 
covers  the  range  of  temperature  of  80  £  T  <  1160°F,  while 
the  set  of  coefficients  No.  2  covers  the  range  of  tempera¬ 
ture  of  1160  £  T  £  2240°F.  The  temperature  of  1160°F  was 
chosen  in  a  way  that  the  standard  deviation  given  by  the 
model  is  almost  the  same  when  either  set  is  used  over  its 
range  as  may  be  seen  from  Table  2. 

It  should  be  noted  that  the  pseudoreduced  temperature 

term  appearing  in  Equation  32  is  defined  as  —I—.  The  pseudo- 

P  c 

critical  temperature,  T  ,  should  be  calculated  using  the 

P  c 

molal  average  not  Equations  66  and  69.  However  if  the  gas 
composition  is  not  available  Equations  66  and  69  may  be  used 
but  higher  errors  should  be  expected. 

Both  sets  of  coefficients  No.  1  and  2  were  obtained 
by  employing  the  Thinh,  et  al  .  equations  (141)  in  obtaining 
the  data  points.  These  sets  of  coefficients  could  be  used 
with  Equations  32,  71,  and  72  for  natural  gases  operating  at 
high  temperatures . 

Equations  33,  73,  74  and  75  were  developed  to  take 
into  account  the  effect  of  presence  of  N2»  C02>  and  H2S. 

These  correction  equations  give  good  results  as  was  indicated 
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earlier  in  Preliminary  Model  Testing;  however,  their  use  is 

subject  to  the  following  restrictions. 

i .  32  <  T  <  300 °F 

ii.  0.0  £  £  0.40  mole  fraction 

iii.  0.0  £  Yc()2  1  0.20  mole  fraction 

iv.  0.0  <  Y m  <  0.25  mole  fraction 

_  11 2  - 

v.  They  may  only  be  applied  when  the  set  of  coefficients 
No.  0  is  used  in  Equations,  32,  71,  and  72. 

2 .  Heat  Capacity  Departure: 

The  proposed  heat  capacity  departure  model  as  expressed 
by  Equation  37  was  compared  with  the  Weiss  and  Joffe  (142) 
generalized  correlation  as  shown  in  Table  5.  When  the  various 
(Cp-C°)  generalized  correlations  (143,144)  were  compared  with 
the  Weiss  and  Joffe  correlation,  the  proposed  model  gave  the 
least  deviation  as  shown  in  Table  5.  This  is  probably  due  to 
the  fact  that  both  correlations  use  the  B-W-R  equation  of  state. 

The  Weiss  and  Joffe  correlation  was  developed  for  pure 
hydrocarbons ,  then  used  to  estimate  the  heat  capacity  depar¬ 
ture  of  oxygen  and  nitrogen.  It  gave  acceptable  results  and 
as  a  result  was  considered  as  a  generalized  correlation 
(145).  In  a  similar  manner  the  proposed  heat  capacity  depar¬ 
ture  model  was  developed  from  an  equation  of  state  which 
describes  the  volumetric  behavior  of  sweet  natural  gases. 
Equation  23,  then  was  tested  by  comparing  with  the  experi¬ 
mental  data  of  Workman  (146),  Krase  and  Mackey  (147)  for 
oxygen  and  nitrogen.  This  model  gives  satisfactory  agreement 
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with  experimental  data  as  shown  in  Tables  3  and  4.  This 
check  was  performed  in  order  to  extend  the  applicability  of 
the  proposed  model  to  non-hy drocarbon  gases  and  thus  to 
sour  natural  gases. 

The  various  comparisons  and  tests  performed  on  the 
proposed  model  suggest  that  it  could  be  used  as  a  genera¬ 
lized  correlation  for  heat  capacity  departure.  This  is  in 
spite  of  the  fact  that  the  coefficients  of  Equation  37  were 
originally  developed  using  only  sweet  natural  gas  compress¬ 
ibility  factor  data. 

It  should  be  mentioned  that  although  errors  in  (C  -  C°) 

P  P 

may  at  times  be  large,  the  effect  on  overall  accuracy  of 

is  usually  much  smaller.  Since  the  value  of  C  is  the  quan- 

C  -  C°  p 

tity  desired,  it  is  the  error  in  (—^7 - which  is  of  import- 

P 

ance  not  the  error  in  (C  -  C°)  (148)  as  shown  in  Table  6. 

p  p 

The  proposed  (C  -  C°)  generalized  correlation  has  an 

r  r 

advantage  over  the  others  in  that  it  is  adaptable  to  computer 
usage . 

The  Hydraulic  Friction  Factor 

The  friction  factor  in  the  rough  pipe  region  is  repre¬ 
sented  by  the  Von  Karman  formula.  Colebrook  (149)  derived  an 
approximation  for  the  friction  factor  in  this  region  when 
the  absolute  roughness  of  pipe  is  known.  The  friction  factor 
equations  suggested  by  Cullender  and  Smith  (150)  are  but  the 
Colebrook  approximation  when  the  absolute  roughness  is  0.0006 


inch. 
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The  assumption  of  fully  developed  turbulent  flow  of 

gas  through  tubing  may  be  in  error.  As  for  example,  for  the 

12  test-points  used  in  testing  main  models  AA  and  BB,  flow 

was  in  the  transition  zone  as  the  comparison  between  Reynolds 

number  (R  )  and  the  boundary  Reynolds  number  (R  ) 
e  eb 

in  Table  7  may  indicate.  As  a  result  the  Cullender  and 

Smith  equations,  which  employ  the  assumption  of  fully  turbu¬ 
lent  flow,  may  give  errors  in  friction  factors  estimated 
ranging  from  -2.17  to  -14.51  per  cent  as  shown  in  Table  7. 
This  order  of  error  may  produce  errors  of  -7.9  to  -0.2  psi 
in  the  calculated  flowing  bottom-hole  pressure  as  the  com¬ 
parison  between  methods  3  and  4  in  Table  12  may  indicate. 

Overall  Heat-Transfer  Coefficient 

In  the  development  of  the  U-j  model  it  was  assumed  that 
the  various  solid  layers  were  tightly  fitted  together  with 
no  intervening  "air  spaces".  If  the  layers  touch  each  other 
only  at  several  points,  the  resistance  to  heat  flow  will 
increase  remarkably  (151).  This  suggests  that  the  thermal 
resistance  of  the  surroundings  in  gas  well  systems  will  have 
a  comparatively  large  value  due  to  the  presence  of  the  annu¬ 
lar  space. 

In  shallow  gas  wells  where  there  exists  one  set  of 
casing  only,  the  composite  cylindrical  wall  consists  of  tub¬ 
ing,  annular  space,  casing,  and  cement.  In  deep  wells  the 
geometry  of  the  composite  wall  is  so  complicated  due  to 
several  factors  such  as:  the  presence  of  more  than  one  set 
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of  casing,  eccentricity  of  both  tubing  and  casing  sets, 
lack  of  cement  at  some  sections,  inevitable  reduction  in 
well  diameter  as  drilling  advances,  caving  at  some  shale 
sections  and  enlargement  at  others,  presence  of  fluid  such 
as  water  in  contact  with  either  casing  or  cement  at  some 
unknown  sections  along  the  well  course,  corrosion,  etc... 

Due  to  such  complexity  and  lack  of  information,  it  is  rea¬ 
sonable  to  assume  that  the  system  has  an  average  thermal 
resistance,  R*,  which  will  be  determined  by  a  trial-and- 
error  procedure.  Furthermore,  this  complexity  and  lack  of 
information  show  it  is  not  possible  to  employ  model  U-j  in 
gas  well  systems. 

In  model  ,  the  term  pj—  is  usually  small,  costly 
to  compute,  and  of  the  order  of  1  to  5  per  cent  of  R*. 

This  suggests  that  model  is  preferable  to  since  it 
requires  considerably  less  computing  time  and  does  not 
affect  the  final  results. 

1  .  Stanton  Number: 

Comparison  of  the  values  of  the  Stanton  number  esti¬ 
mated  using  the  various  models  indicates  that  there  is  lit¬ 
tle  difference  for  low  Prandtl  Number  fluids  like  air  (152). 
However,  the  best  expression  that  reproduces  experimental 
data  is  that  of  Gowen  and  Smith  (153).  The  reliability  of 
this  model  in  predicting  the  Stanton  numer  lies  in  an 
implicit  assumption  that  the  friction  factor  and  hence  the 
equivalent  sand  relative  roughness  (— )  completely  charac¬ 
terizes  the  effect  of  roughness  on  temperature  profile. 
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2 .  Prandtl  Number  of  Natural  Gases  : 

The  Prandtl  number,  which  is  necessary  to  calculate 
the  Stanton  number,  is  defined  as  in  Equation  50.  For  a 
gas  mixture,  it  is  recommended  that  the  Prandtl  number  be 
calculated  using  Equation  50;  however,  the  various  terms 
appearing  in  the  definition  should  be  estimated  separately 
(154). 

Thermal  Conductivity  of  Gas  Mixtures 

Equation  56  applies  to  paraffins  giving  an  average 
deviation  of  2.6  per  cent  from  the  experimental  data  (155), 
while  Equation  57  applies  to  non-hydrocarbons  giving  an  aver¬ 
age  deviation  of  2.2  per  cent  from  the  experimental  data 
(156).  The  Sutherland  constant  estimated  from  boiling  point 
temperature  as  expressed  by  Equation  54  may  not  be  reliable; 
however,  an  error  of  20  per  cent  in  the  Sutherland  constant 
produces  only  an  error  of  1  per  cent  in  the  estimated  gas 
mixture  thermal  conductivity  at  atmospheric  pressure  (157). 
Hence  the  simplification  introduced  by  Equation  54  (due  to 
necessity)  is  justified. 

Equations  59,  60,  and  61  for  the  estimation  of  the 
thermal  conductivity  departure,  are  usually  applicable  to 
pure  non-polar  gases.  However,  Reid  and  Sherwood  (158)  have 
recommended  the  use  of  these  equations  for  gas  mixtures. 

Their  recommendation  was  based  on  an  assumption  that  a  gas 
mixture  can  be  considered  as  a  hypothetical  pure  gas  with 
pseudocri  ti  cal  properties.  This  suggestion  has  been  accepted 
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because  few  experimental  data  for  thermal  conductivity  of 
high-pressure  gas  mixtures  are  available. 

The  gas  mixture  thermal  conductivity  model  is  expected 
to  give  errors  of  less  than  5  per  cent  for  simple  non-polar 
gas  mixtures  for  Tr  <  1.3.  The  same  model  may  be  used  for 
polar  and  non-pol  ar polyatomi c  gas  mixtures,  but  larger  errors 
may  occur  (159). 

Viscosity  of  Pure  and  Natural  Gases 

Equation  63  was  originally  developed  for  monatomic 
gases.  However,  it  is  found  to  be  remarkably  good  for  poly¬ 
atomic  gases  as  well  (160). 

Equation  64  is  an  approximation  for  the  collision 
integral  function,  ft  ,  for  non-polar  gases.  It  deviates  1 
to  2  per  cent  from  the  tabulated  function  (161).  It  may  also 
be  applied  to  polar  gases;  however,  errors  up  to  2  per  cent 
may  occur  (162).  The  collision  integral  function  is  required 
for  the  estimation  of  pure  gas  viscosities  at  atmospheric 
pressure  as  given  by  Equation  53. 

The  Lee,  et  a  1 .  correlation  (163)  for  calculating  the 
viscosity  of  sweet  natural  gases  gives  errors  ranging  from 
-11.05  to  +16.38  per  cent  from  their  exoerimental  data.  The 
presence  of  N^,  CO^,  and  tends  to  increase  the  viscosity 

of  natural  gases  (164).  As  a  result  the  Lee,  et  al.  correla¬ 
tion  may  give  errors  larger  or  smaller  than  those  mentioned 
when  applied  to  sour  natural  gases. 
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Flowing  Bottom-Hole  Pressure 

As  may  be  seen  from  Tables  12  and  13,  method  1  always 
underestimates  the  bottom-hole  pressure.  Methods  2  and  4 
underestimate  11  of  12  test-points,  while  method  3  underesti¬ 
mates  9  of  12  test-points.  It  can  also  be  seen  that  method 
3  deviates  the  least  from  the  measured  bottom-hole  pressure 
for  10  of  12  test-points. 

Table  14  indicates  that  methods  2,  3,  and  4  estimate 
pressure  drops  for  test-points  ,  Ag,  A^ ,  B-j  ,  C,  and  D  with 
accuracy  comparable  to  those  of  the  measured  values.  The 
various  methods  estimate  flowing  bottom-hole  pressures  in 
excellent  agreement  with  the  measured  values  for  test-points 
A^,  C,  and  D  as  shown  in  Table  14.  This  is  probably  due  to 
using  accurate  surface  pressure  gauges  (Dead  Weight  Gauges) 
in  test-points  C  and  D. 

The  analysis  of  Table  14  exhibits  what  may  be  a  syste¬ 
matic  error  in  the  calculated  bottom-hole  pressure  as  the 
flow  rate  increases  as  is  shown  in  Figure  3.  This  behavior 
may  explain  the  excellent  agreement  of  the  various  methods 
with  test-point  A ^ . 

Figure  4  shows  examples  of  estimated  (methods  1  and  2) 
and  assumed  (methods  3  and  4)  temperature  profiles.  It  shows 
that  the  fluid  temperature  at  every  depth  is  decreasing  in 
the  following  order:  methods  1,  2,  and  either  3  or  4.  It 
can  also  be  seen  that  the  fluid  temperature  is  always  great¬ 
er  than  or  equal  to  the  earth  temperature.  This  suggests 
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FLOW  RATE  Qq  (mmscf/day) 


FIGURE  3 


DEVIATION  OF  PRESSURE  DROPS  CALCULATED 
BY  METHOD  2  FROM  MEASURED  VALUES  FOR 
ALL  TEST  POINTS 


TEMPERATURE  (°F) 
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DEPTH  (thousand  feet) 


FIGURE  4.  FLUID  TEMPERATURE  PROFILES  CALCULATED 
BY  METHODS  1,  2,  3,  AND  4  FOR  TEST- 
POINT  C 
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that  in  these  examples  heat  is  being  transferred  from  the 
flowing  fluid  to  the  surroundings.  Furthermore,  it  sug¬ 
gests  that  methods  3  and  4  are  equivalent  to  a  hypothetical 
case  where  heat-transfer  is  such  as  to  bring  the  fluid 
temperature  profile  to  a  straight  line. 

The  analysis  of  Table  12  indicates  that  for  any  test- 
point,  method  1  predicts  the  lowest  flowing  bottom-hole 
pressure  followed  by  methods  2,  4,  and  finally  method  3. 

This  behavior  is  consistent  in  all  test-points  for  methods  1, 
2,  and  3.  However  at  three  test-points,  namely  Ag,  Ag  ,  and 
D,  methods  2  and  4  switch  positions.  This  may  be  due  to  the 
use  of  the  Cullender  and  Smith  friction  factor  equations 
(which  deviate  from  the  Von  Karman  rough  pipe  law  by  -0.04%, 
-0.04%,  -0.41%  for  test-points  Ag ,  Ag,  and  D  respectively) 
in  method  4. 

Reasons  for  Deviations  of  Calculated  from  Measured 

Bottom-Hole  Pressures' 

This  discrepancy  between  measured  flowing  bottom-hole 
pressures  and  those  calculated  by  method  2  may  be  attributed 
to  a  combination  of  the  following  factors: 

1  .  The  Existence  of  Two-Phase  Flow: 

The  presence  of  condensate  in  the  flowing  gas  has 
been  dealt  with  according  to  Standing  (165)  when  the  gas/ 
condensate  ratio  is  greater  than  40,000  SCF/bbl.  He  suggests 
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that  the  gravity  of  the  flowing  stream  (Gf)  be  used  in¬ 
stead  of  that  of  dry  gas  ( G )  according  to: 


G  +  4.591  P-C-°nd- 

af  =  ,  '  'TT2'3' - 

i? 


(71) 


where  R  =  Gas/Condensate  Ratio,  SCF/bbl. 

3 

and  pcond  =  density  of  condensate,  gm/cm  . 

The  pseudocri ti cal  pressure  and  temperature  of  the  flowing 
stream  may  then  be  estimated  using  G^. 

This  method  of  dealing  with  two-phase  flow  does  not 
account  for  the  effect  of  the  liquid  phase  presence  on  the 
isobaric  heat  capacity,  viscosity,  thermal  conductivity, 
and  friction  factor. 


2 .  The  Presence  of  Water  Vapor  in  the  Flowing  Fluid: 

The  presence  of  water  vapor  in  the  flowing  stream  has  not 
been  accounted  for.  Although  the  water/gas  ratio  is  small  (1.5 
bbl/MMSCF)  except  for  test-point  C  (5.2  bbl/MMSCF),  water  is 
highly  polar  and  as  a  result  the  behavior  of  the  water-gas  mix¬ 
ture  system  may  deviate  from  the  law  of  corresponding  states. 
This  behaviour  may  not  be  estimated  accurately  by  the  correla¬ 
tions  employed  and  as  a  result  errors  in  the  estimated  pres¬ 
sure  may  be  introduced,  However,  it  is  surprising  to  find  that 
test-point  C  is  in  very  good  agreement  for  the  various  methods. 

3.  Well  Bore  Deviation  from  Vertical: 


Due  to  lack  of  information  to  the  contrary  all  calcula- 
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tions  of  flowing  bottom-hole  pressure  were  made  assuming 
that  the  tubing  was  situated  in  the  vertical  position.  How¬ 
ever,  in  practice  well  bore  and  therefore  tubing  segments 
may  deviate  from  the  vertical  as  is  shown  in  Figure  5.  Such 
deviations  may  result  in  a  decrease  in  pressure  drop,  in 
which  case  the  calculated  flowing  bottom-hole  pressure  may 
show  positive  deviations. 

4 .  Unreliable  Estimates  of  the  Viscosity  of  Sour  Gases: 

The  presence  of  appreciable  amounts  of  CO^,  N£,  and 
H^S  tends  to  increase  the  viscosity  of  hydrocarbon  mixtures 
(167).  The  Lee,  et  al .  (168)  viscosity  correlation,  there¬ 
fore,  may  underestimate  or  overestimate  the  viscosity  of 
sour  natural  gases  depending  on  the  sign  and  magnitude  of 
errors  produced  by  this  correlation  when  applied  to  hydro¬ 
carbon  mixtures.  Consequently  Reynolds  number  and  hence 
the  friction  factor  may  be  in  error.  Such  errors  may  be 
significant  especially  at  relatively  low  Reynolds  number. 


5 .  Uncertainty  of  the  Length  of  Flow  String: 

The  measured  flowing  bottom-hole  pressures  were  reported 
at  the  mi d-poi nt-of-perforati ons  (MPP).  It  is  not  known 
from  the  available  data  whether  the  measuring  device  was  loc¬ 
ated  at  that  point  or  whether  in  fact  the  tubing  extended  to 
the  mid-point  of  perforations.  If  the  tubing  does  not  extend 
to  the  mi d-poi nt-of-perforati ons  ,  method  2  tends  to  under- 
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FIGURE  5.  PLAN  VIEW  OF  14  VERTICAL  WELLS 
DRILLED  TO  6000  FEET  (166) 
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estimate  the  flowing  bottom-hole  pressure. 

6 .  Variation  of  the  Overall  Heat-Transfer  Coefficient: 

Models  U 2  and  (due  to  necessity)  were  employed  as 
a  means  for  estimating  the  overall  heat-transfer  coefficient 
for  complex  geometries  such  as  those  encountered  in  gas  wells. 
However,  these  approximations  may  not  be  reasonable  from  a 
physical  standpoint.  Consequently  the  heat- trans fer  between 
the  flowing  fluid  and  the  surroundings  may  be  in  error  which 
may  affect  the  calculated  flowing  bottom-hole  pressure. 

Main  models  AA  and  BB  employ  models  and  respect¬ 
ively.  It  has  been  found  that  these  main  models  yield  cal¬ 
culated  pressure  and  temperature  profiles  which  are  essenti¬ 
ally  the  same;  however,  the  cost  of  computation  may  be 
reduced  by  about  a  factor  of  10  when  using  main  model  BB 
instead  of  AA  as  may  be  seen  from  Table  15.  For  example  in 
six  test-points  the  range  of  cost  for  the  calculations,  when 
the  computer  programs  were  stored  on  disks,  was: 

Main  model  AA  $3.18  -  $10.83 

Main  model  BB  $0.44  -  $  0.75 

7 .  Uncertainty  of  Absolute  Roughness  of  Tubing: 

This  is  probably  one  of  the  most  important  sources  of 
error.  The  absolute  roughness  of  tubing  was  taken  as  0.0006 
inch.  This  value  is  that  shown  by  Cullender  and  Binckley  (169) 

for  new  tubing.  This  value  would  in  all  likelihood  be  too  low 
for  old  pipe  as  indicated  by  Colebrook  (170).  He  reported 
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that  an  increase  in  roughness  may  cause  20  to  30  per  cent 
reduction  in  the  carrying  capacity  of  pipes  while  the  reduc¬ 
tion  in  the  area  perpendicular  to  flow  due  to  such  increase 
in  roughness  may  cause  a  reduction  of  only  2  to  3  per  cent 
of  the  pipe  carrying  capacity.  This  implies  that  the  same 
increase  in  roughness  may  cause  an  increase  of  25  to  43  per 
cent  of  pressure  loss  due  to  friction  if  the  pipe  carrying 
capacity  were  kept  at  its  initial  value.  Furthermore,  Ippen 
(171)  has  reported  a  practical  example  where  the  absolute 
roughness  was  doubled  within  three  years  and  the  associated 
friction  factor  was  increased  by  20  per  cent. 

This  suggests  that  flowing  bottom-hole  pressure  calcu¬ 
lated  using  the  value  of  roughness  of  0.0006  inch  will  yield 
values  which  are  usually  low.  The  12  test-points  presented 
here  tend  to  support  this. 
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CONCLUSIONS 


As  a  result  of  this  investigation,  the  following  con¬ 
clusions  can  be  made : 

1.  The  assumption  of  adiabatic  flow  in  gas  wells  is 
unrealistic  and  yields  the  lower  bound  estimate  for 
flowing  bottom-hole  pressure. 

2.  The  apparent  failure  of  main  models  AA  and  BB  may  be 
attributed  to  unreliable  data  and  errors  in  the  differ¬ 
ent  terms  such  as  viscosity  and  absolute  roughness, 
etc.  .  . 

3.  The  complete  solution  of  the  Force-Momentum-Balance 
differential  equation,  as  suggested  in  this  study, 
yields  the  closest  estimates  to  the  measured  flowing 
bottom-hole  pressures. 

4.  The  assumption  of  absolute  roughness  of  0.0006  inch 
for  old  and  new  tubing  may  be  in  error.  In  all  like¬ 
lihood  this  value  is  too  small. 
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RECOMMENDATIONS 


The  following  recommendations  are  made: 

1.  The  Lee,  et  al.  (172)  viscosity  correlation  should  be 
replaced  by  a  more  accurate  one  capable  of  represent¬ 
ing  the  viscosity  of  sour  as  well  as  sweet  natural 
gases.  In  this  regard  it  is  recommended  that  the  Carr, 
et  al.  (173),  Dean  and  Stiel  (174)  viscosity  correla¬ 
tions  be  examined. 

2.  The  parameters  in  the  equations  giving  the  pseudocri- 
ticals  of  natural  gases  as  a  function  of  gas  gravity 
should  be  re-estimated  by  using  data  for  real  natural 
gases  rather  than  fitting  the  Brown,  et  al .  (175) 
correlation. 

3.  The  absolute  roughness  of  tubing  used  in  gas  wells 
should  be  re-estimated.  This  can  be  done  by  employing 
accurate  surface  and  bottom-hole  measurements  for 
large  number  of  dry  gas  wells  and  applying  main  model 

A  or  B.  The  values  of  absolute  roughness  estimated 
may  then  be  correlated  to  age  of  tubing,  the  degree  of 
sourness  of  flowing  gas,  etc... 
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NOMENCLATURE 


A,  B 

=  functions  of  snecified  arguments  or 
constants  depending  uoon  usage. 

B* 

=  a  function  of  specified  argument. 

c 

P 

C* 

=  isobaric  heat  capacity,  Btu/lbm.°R. 

=  a  function  of  specified  argument. 

o 

Cp 

=  molal  ideal  isobaric  heat  capacity, 

Btu/1 b  mol e . °  R . 

o 

Cpm 

m 

=  molal  ideal  isobaric  heat  capacity  of 
gas  mixtures,  Btu/lb  mole.°R. 

O 

Cp 

mH .  C  . 

=  molal  ideal  isobaric  heat  capacity  of  pure 
hydrocarbon  mixtures,  Btu/lb  mole.°R. 

o  > 
o 

=  molal  specific  heat  at  infinite  volume, 
Btu/lb  mol e . °  R . 

Cp 

=  molal  isobaric  heat  capacity,  Btu/lb 
mol e . °  R . 

Cp 
r  m 

=  molal  isobaric  heat  capacity  of  gas 
mixtures,  Btu/lb  mole.°R. 

C 

V 

=  molal  constant-volume  specific  heat, 

Btu/lb  mol e . °  R . 

d 

=  inner  pipe  diameter,  inches. 

D 

=  inner  pipe  diameter,  feet. 

D* 

=  a  function  of  specified  argument. 

Do 

=  inner  pipe  diameter,  feet. 

exp 

=  exponential. 

E* 

=  a  function  of  specified  argument. 

ep 

=  estimated  truncation  error  in  calculated 
pressure,  1 b / i n  2 . 

et 

=  estimated  truncation  error  in  calculated 
temperature,  °F. 
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bi 

. 

_ 

. 

' 
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f 

f  1  »  f  2  >  fen 


FH2S’  fn2 

g 


9C 

G 


G 


G 


f 


h 


h 


o 


J 

k 

k 

k° 


Moody's  friction  factor,  Dimensionless, 
functions  of. 


functions  of  specified  arguments. 


local  gravitational  acceleration, 

32.174  f t/sec^ . 

p 

conversion  factor,  32.174  1 bm. f t/1 . sec  . 

molal  flow  rate  per  unit  area, 
lb  mol e/f t^ . sec . 

gas  gravity,  (air  =  1 ) . 

gravity  of  flowing  stream, 

(air  =  1 ) . 

enthalpy,  Btu/lbm. 

gas-film  heat-transfer  coefficient, 

Btu/ft^ . sec . °R. 

trial  step  size  at  kth  interval,  feet. 

conversion  factor,  778.16  lb^.ft/Btu. 

upper  limit  integer. 

thermal  conductivity,  cal /cm. sec . K. 

thermal  conductivity  at  atmospheric 
pressure,  cal /cm. sec . K. 

thermal  conductivity  of  gas  mixtures,  at 
atmospheric  pressure,  cal/cm.sec.K. 

thermal  conductivity  of  gas  mixtures, 
cal /cm. sec . K. 


■ 
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K 

In 


log 

L 

L* 


M 


W 


n 

P 


pPc 


P  c 


H.C. 


q 

Q. 


R 


functions  of  specified  arguments. 


a  function  of  specified  argument. 

natural  logarithm. 

logarithm,  base  10. 

length  variable,  feet. 

length  of  tubing,  feet. 

molecular  weight  or  average  molecular 
weight  depending  upon  usage,  lbm/lb  mole. 

upper  limit  integer. 

2 

pressure  variable,  lb^/ft  absolute. 

critical  pressure  or  pseudocritical 
pressure  depending  upon  context, 

Ib^/in^  absolute. 

2 

pseudocritical  pressure,  lb^/in  absolute 

pseudocritical  pressure  of  pure  hydro¬ 
carbon  mixtures,  Ib^/in^  absolute. 

reduced  or  pseudocritical  pressure  depend 
ing  upon  context,  Dimensionless. 

Prandtl  Number,  Dimensionless. 

heat  flux,  Btu/sec . 

flow  rate  at  standard  conditions, 
MMSCF/Da.y . 

universal  gas  constant, 

1545  lbf.  ft/lb  mole  °R. 
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R 

R* 


R** 


R 


e 


R 


e 


b 


S 


S 

T 


T 


b 


T 


c 


P  c  H  .  C  . 


T 


r 


T 


s  r 


u 


<u> 

U 


V 

w 


=  Gas/Condensate  Ratio,  SCF/bbl. 

=  average  thermal  resistance  of  the  surround¬ 
ing  as  defined  by  Equation  13, 
f t^  .  sec  . 0  F/Btu  . 

=  average  overall  heat-transfer  resistance, 
f t2 . sec . 0  F/ Btu  . 

=  Reynolds  Number,  Dimensionless. 

=  Reynolds  Number  at  the  boundary  between 
the  transition  and  rough  pipe  regions, 
Dimensionless. 

=  Sutherland  Constant,  K. 

=  sin  6,  dimensionless. 

=  temperature  variable,  °R. 

=  normal  boiling  point  temperature,  K. 

=  critical  or  pseudocritical  temperature 
depending  upon  context,  °R. 

=  pseudocritical  temperature,  °R. 

=  pseudocritical  temperature  of  pure  hydro¬ 
carbon  mi xtures  ,  °  R  . 

=  reduced  or  pseudoreduced  temperature. 
Dimensionless. 

=  temperature  of  surroundi nas  ,  °R. 

Ts 

=  j—  ,  Dimensionless, 
c 

=  velocity  of  flow,  ft/sec. 

=  average  velocity  of  flow,  ft/sec. 

=  overall  heat-transfer  coefficient, 

Btu/ f t^ . sec . 0  F . 

=  molal  volume  variable,  ft^/lb  mole. 

=  shaft  work  done  by  gas,  lb^.ft. 
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X,  Y 
Y 

yco2 

yh2s 


z 


Z 


z 


c 


=  functions  of  specified  arguments. 

=  mole  fraction.  Dimensionless. 

=  mole  fraction  of  carbon  dioxide. 

=  mole  fraction  of  hydrogen  sulphide. 

=  mole  fraction  of  nitrogen. 

=  height  above  datum,  feet. 

=  compressibility  factor.  Dimensionless. 

=  critical  compressibility  factor. 
Dimensionless. 


Greek 


a 

Y 

6 

6 1  >  6  2 
AP 

e 


e 


* 


£  3 


0 

y 


=  inverse  of  the  momentum  correction  factor. 

=  a  parameter  defined  by  Equation  58. 

=  a  change  in  the  value  of  a  variable. 

=  tolerances  used  as  conversion  criteria. 

2 

=  pressure  drop,  lb^/in 
=  absolute  roughness,  inches. 

=  /~t~  •  R  •  §  =  Shear  Reynolds  number  x 

Bed 

relative  roughness. 

=  maximum  energy  of  attraction  for  the 
Lennard-Jones  potential,  erg. 

=  critical  temperature  correction  as  defined 
by  Equation  24,  °F. 

=  equivalent  sand  roughness,  inches. 

=  dip  angle  of  pipe,  degrees. 

=  viscosity,  c.p. 
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V 

y 

p 

p 

p 

p 


m 

c 

cond 

r 


a 


I 

X 


=  viscosity  at  one  atmosphere,  c.p. 

=  viscosity  of  gas  mixtures,  c.p. 

3 

=  density  variable,  Ibm/ft  . 

3 

=  critical  density,  lbm/ft  . 

3 

=  density  of  condensate,  gm/cm  . 

=  reduced  density,  Dimensionless. 

=  collision  diameter  for  the  Lennard- Jones 
Potent i al ,  & . 

=  summation  operator. 

=  a  function  of  specified  argument  as 
defined  by  Equation  49. 

=  collision  integral. 
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APPENDIX  A 


DEVELOPMENT  OF  A  MATHEMATICAL  MODEL  FOR  STEADY 
STATE  COMPRESSIBLE  FLOW  OF  GAS  IN  PIPES 


For  the  steady  state  flow  of  gas  through  a  differen¬ 
tial  length  of  pipe  of  known  dimensions  and  position  in 
space,  the  Total  Energy  Balance  Equation  and  the  Force 
Momentum  Balance  Equation,  based  on  the  fluid  element  con¬ 
tained  in  a  differential  length  (dL)  may  be  written. 

The  Total  Energy  Balance  Equation: 


6q 


dW 

J 


tiD2  G 


dh  + 


<u>d<u>  + 
agcJ 


JL  dz 


(A-l) 


The  Force  Momentum  Balance  Equation: 


(BT-JdP 


<u>d<u> 


(A-2) 


where  <5q 


W 

G 

h 

<u> 

P 


the  transferred  heat  from  the  surroundings  to 
the  fluid  element,  Btu. 

shaft  work ,  lb^.  ft. 

2 

molal  flow  rate  per  unit  area,  lb  mole/sec.  ft. 

enthalpy,  Btu/lbm. 

mean  velocity  of  fluid,  ft/sec. 

2 

pressure,  lb^/ft  absolute. 
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3 

V  =  molal  volume,  ft  /lb  mole, 
f  =  Moody's  hydraulic  friction  factor, 
z  =  elevation  above  datum,  ft. 

D  =  I . D.  of  pi pe ,  ft. 

Mw  =  the  molecular  weight  of  gas,  Ibm/lb  mole, 
a  =  the  inverse  of  the  momentum  correction  factor 


J  =  conversion  factor. 

=  778.16  lbf.ft/Btu 

g  =  local  gravitational  acceleration. 

=  32.174  ft/sec. ^  (assumed). 

g  =  conversion  factor, 
c 

Ibm.ft , 

=  32.1  74  — - - * 

1  b ^ .  se cc 

Equations  ( A - 1 )  and  ( A- 2 )  were  rearranged  into  forms 
appropriate  to  numerical  solution.  This  task  was  accom¬ 
plished  as  follows: 

By  definition: 


Thus 


<u> 
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<u>  =  GV 


(A-  3) 


and 
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uaavtt.^dr  = 

.  no  id g'l 9 T 90D6  Fsrio  rdsd  rvfiTg  FbdoF  =  g 

.  (b9fflU226)  \D92\dd  .$£  a 

. TOdoed  oe  r  ^sv.ioo  a 

r— - m.se  - 


annol  odm  b9gn6't"<69't  s^ew  (S-A)  bne  (F-a)  enordsup? 

.  n  o  r  J  u  f  o  2  Fsoribrnuti  od  sd  sq  j  tqqB 

: awoT  Fod  26  b9dar  Fq 


:  '  Jf  rti  tab  y,8 


<u> 


a 


eurtT 


va  -  <u> 


bns 


Vb  3  =  <u;  j 
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The  dip  angle  of  the  pipe  is  0,  where 


sin  0  =  s  = 


dz 

dL 


( A  -  5 ) 


The  equation  of  state  is: 


PV  =  ZRT 


( A  -  6 ) 


From  which 


dV  =  -y-  dT  +  ^  dZ  -  dP 


(A- 7 ) 


and 


,d\_\  ZR  RT  ,3Zv 

V9T'P  P  P  ^T;P 


(A-8) 


where  Z  =  /<jh(T  ,Pr).  Thus 


dz  -  (gT  )p  dT  +  (gp  )T  dp 
r  r  r  r 


or 


dz  =  (|f)p  dT  +  (II) T  dP 


SP'T 


( A  -  9 ) 


The  gas  does  no  shaft  work;  that  is. 


dW  =  0 


(A- 10) 


9T9I1W  ,9  z\  aq  |  artt  jfgr.fi  qrb  9riT 


(2-A) 


:2'r  eteta  to  norjfiupa  artT 


TflS  =  VS 


rtsfriw  men  -i 


+  Tb 


rcs 


,  ,  Tfj  ,  SS  _  ,  V_.  i 

SlT6'  S 


sriw 


*t  o 


c-'  T(ff)  *  Tb  „(!!)  •  Tb 


qvTG 


z\  JerlJ  ;;>how  :}tfi(l2  o.i  29ob  asp  9fiT 


0  =  Wb 
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The  transferred  heat  to  the  fluid  element  from  the  surround¬ 
ings  may  be  given  by 


Sq  =  -uDdL  U(T  -  T) 

3 


(A-ll) 


where  T  =  the  temperature  of  the  surroundings,  °F. 

T  =  the  temperature  of  the  flowing  fluid  element,  °F 

? 

U  =  the  overall  heat- trans fer  coefficient,  Btu/ft  . 

°  F . 


It  is  known  that  h  =  fen( T,P).  Thus 


dh  =  ($)p*dT  +  (|^)T-dP 


dh  =  Cp-dT  +  (|^)T-dP 


(A- 1 2) 


But  ^3 P  ^  T 


V  -  T(|^) 


3T '  P 


JM 


since  h  is  in  Btu/lbm  units.  ( A- 1 3 ) 


w 


Equations  (A-6),  (A-7),  ( A - 1 2 )  ,  and  ( A- 1 3 )  give 


RT1 


dh  -  cp-dT  -  JpM  v3T,p 

w 


(Ir)D-dP 


(A- 1 4) 


where  cp  is  in  Btu/lbm.  °R  units 


-bmjo't'tuz  9fi.j  mo^t  Jrisir-  f  ■  b  ot  £.v  eTist  2(16*1$  sriT 

yc  rovig  91  \^6m  sgnf 


T )  U  JbQir  =  pj 


.1°  t  egnrbnuo't^ue  srit  to  s'luis't 9 qcr> JJ  s  i  «  T  snariw 
,i!9m9r9  aruTt  gnrwort  b\M  to  s^ut  6T9q>i9J  sri  i  = 

.  Jt\uJ8  ttn9'roitt900  *\9t  ane^t-tsart  lU^syo  srIJ  *  U 

aurlT  .  ( ^ «  T ) =  rl  tsdt  nwoml  ar  JI 


qb-T(f )  +  Tb-^'li)  -  ,1b 


.jJrnu  mdf\uJ8  nr  at  rl  9onra 

w 

.2trnu  .mdr\ut8  nf  z\  qD  s^eriw 


T  ‘  96^ 


i  u8 


,  (  <  - '  2fl0  Jsup3 
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<u>d<u>  G  VdV 


ag 


ag, 


by  virtue  of  Equations  ( A - 3 )  and  ( A- 4 ) . 


G2Z2R2T2  / dT  dZ 

- ? —  ‘  T*  T 

“9cP 


dPx  with  the  aid  of  /„  lc\ 
T}  Equation  (A-7).  (A'15) 


f<u>2  dL  =  f G2Z2R2T2  .  ..  with  the  aid  of 

29CD  2g  DP2  Equation  ( A- 3 ) 

c 


( A- 1 6 ) 


V  I D  _  ZRT  ,D  with  the  aid  of 

M  K  PM  Equation  (A-6) 

w  w 


( A- 1 7 ) 


q  .  q  _  byvirtueof 

g  dz  ~  g  dL  Equation  ( A - 5 ) 


(A- 18) 


Combining  Equations  (A-l),  (A-9),  ( A- 1 0 ) ,  ( A- 1 1 ) , 
(A-14),  (A-15),  and  (A-18)  and  rearranging  yield 


dP 

dL 


G2ZR2T2  (dl. 


RT1 


agcJP 


2  V9P'T  J  PM 


w 


(1Z) 

l9T'P 


G2Z2R2T2 


agc  J  P 


dT 

dL 


G2Z2R2T  G2Z2R2T2  /3Zv 
CP  nn2  ln2  ^T'P 


«gcJP' 


agcJP‘ 


4u<Ts  -  T)/Hw  -  f-  DGS 

c 

DGJ 


Combining  Equations  (A-2),  (A-9),  (A-15),  (A-16), 


( A- 1 9  ) 
(A- 17) 


and  (A-18)  and  rearranging  yield 


.(£-A)  brts  (S-A)  srtor^fup. 


. ( K-A)  no r leup3 


Jb 


"to  br&  ari*  rttrw 


(8f-A) 


,  *?  suiH/v  :d  JL  =  sb  X 


«(IT-A)  « (or-A)  «(G-A)  « ( r -a)  enordfupa  gnrnrdmoD 

bfsrx,  onf’gnfi't'tBSi  bn&  (8  -A)  bn^  ,  Sf-  '  ,(M-A) 


qvT6/ 


'■  1“ )  - 


1  3  l  o , 


,  ...  SS%S3  T  fl 


«(G-A)  «  ( S ~ A )  anorjsup’-t  jnrnrdmoO 


bie;  pnrgn  b . ■£■  (8f-A)  ins 
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dP 

ai 


G2RTP  ,3Zn 
a  l3P'T 


G2ZRT 

a 


+ 


dT  G2ZRP  ,  G2RTP  ,3Zn 
dL  a  a  ^T;P 


9cP  5g  fG2ZRTP 
ZRTgc  '  2D 


(A-20) 


By  employing  the  definitions  of  Pr  and  and  multiply¬ 
ing  both  sides  of  Equation  ( A - 1 9 )  by  -1,  Equations  ( A - 1 9 ) 
and  (A-20)  may  be  written  as: 


dP 


d L 


G2ZR2T2  ,3Z  >  RTTr  ,3Z  x  ,  G2Z2R2T2 

^77^  JPrMw  C3Tr>Pr  — 


dT, 

dL 


X  -  .  G2Z2R2T2 

Tr  P  ag  JP2T 


G2ZR2T2  ,3Z  x 
777"  3Tr  Pr 


and 


f-  DGS  -  4UJ<Tsr  -  V/Mw 


DGJ 

dPr 

g  p3 

G2RTZP 

G2RTP  , 

dL 

P  M 

aPw 

a  ^ 

r  w 

r 

dTr 

G2ZRPT 

G2RTP 

(  9Z  ) 
l9T/P. 

dL 

aT 

T  — - — - - 

a 

r 

r  r 

gSP3 

fG2ZRTP 

ZRT 

2D 

3P  T 
r  r 


(A-21) 


( A- 2  2 ) 


w 


(OS -A) 


qjqsVt 


?TrS 


/ c  <\nh  ,  irs'-d 


-Y.rqr^rum  bns  T  brtfi  3  "to  ano f flr^sb  9rtJ  r y o f qs.: e3  y8 

( 6 f -/  )  . 2 r  > i  p  I  t  F -  yd  ( e f -A)  n o r;  6 u p 3  f-o  ?->  r  riJod  gnr 

:zb  nsJj-Mw  9d  ysm  (OS-A)  bns 


T, 


qc  gjo  -t  * 


vPv^  ,T  r-;;el;6 


Sb 


q  \  jn  ' 


"IF 


WM\(,T  -  -  *00  i 


o 

.n 

~TB 

ql  tg; 

•,Tb 

-  - -  T 

IB 
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and 

where 


Equations  ( A- 2 1 )  and  (A-22)  may  be  written  as 


dP 


dT 


D*  HT1  +  E*  =  F* 


dP  dT 

A*  +  B*  =  C*  respectively 


A*  = 


9cP 

PrMw 


g2rtzp  ,  g2rtp  ,az  x 
a 


aP. 


■  /««(Tr.Pr.Tc.Pc.Mw.Q0.D) 

_  G2ZRPT  ,  G2RTP  t  9Z  v 

B  '  ~^T7  +  a  (5TV)Pr 
'  /»»(Tr.Pr.Tc.Pc.Mw.Q(>.D) 

„*  qSP3  fG2ZRTP 
C  "  "flT - 2D— 

■  /^(Tr.Pr.Tc.Pc,Hw>Q0,D,(f),Re>5) 

But  Re  =  =  /<3«(Tr,Pr>Q0>D) 


(A-23) 

( A- 2  4 ) 

( A- 2  5 ) 

( A  -  2  6 ) 


PS-A) 


Vfi3  9i  2  S'* 


.  1 


P 

"T'[  i 


(0'oP<wM'd<,-3T',‘1«-,T'  ■'«* 


(a<0P-WM'3<l'3T<iq'  * 


<C-,P.W'  ..  -3'  > 


(C.o0,.S._T)«O\  -C  p  . 


9T9riw 


Ju8 
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C*  =  /on(T„>P„>T^,P,,M,.1>Q„,D,(|),S)  (A-27) 


r ’  r  ’  c  ’  c  7  w’^o 


D*  = 


G2ZR2T2,3Z  >  RTTr  ,3Z  >  G2Z2R2T2 

ag  _JP2  3Pr  '  r  JPrMw  (VPr  ag  JpV 


‘  /«»(Tr.Pr.Tc.Pc.Mw.Q0.D) 


(A-28) 


E*  = 


T 


T  CP  ? 

'r  K  ag  JP^T 


G2Z2R2T2  _  G2ZR2T2  / 9Z  v 

ag  JP2  BTr  Pr 


=  /^(Tr»pr»Tc»Pc»Mw»Qo>D)  (A_29) 


-9-  DG  -5  -  4UJ  (T  -  T  )/M 
g  v  s  r  r  w 

 fc 

DG  J 

=  /^n(Tr,Pr,Tc>Pc,Mw,Q0,5,U)  (A-30) 


Solving  Equations  ( A- 2 3 )  and  (A-24)  simultaneously  for 
dPr  dTr 

nr  and  nr  gives 


C*  •  E* 
A*  •  E* 


B*  •  F* 
B*  •  D* 


( A- 3 1 ) 
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and 


A* •  F *  -  C * •  D * 
A* •  E*  -  B*  •  D* 


( A- 3 2  ) 


dP  dT 

In  order  to  obtain  a  solution,  and  must  be 

finite. 

Equations  (A-31)  and  (A-32)  with  the  aid  of  Equations 
(A-25)  through  (A-30)  may  be  rewritten  as 


dP 


^  -  fon(  T  ,P  ,T  ,P  ,^.Q0.D.|.U.J?) 


(A-33) 


and 


dT 


HIT  =  ^(VPr’Tc>Pc>WD-7T’U’5> 


(A-34) 


For  a  specific  flow  case,  that  is,  when  the  pipe  dimen¬ 
sions  and  position  in  space  (L*,D,s')»  flow  rate  (Q0)  and  gas 
properties  (T  ,P  ,M  )  are  given.  Equations  (A-33)  and  (A-34) 

V  V  W 

reduce  to 


dP 


dL 


/cn(Tr,Pr,(f),U) 


( A- 3  5 ) 


and 


dT 


cTT 


r  


fcn( Tr,Pr,(|),U) 


( A-36 ) 


Equations  (A-35)  and  (A-36)  constitute  the  mathematical 


model  for  a  steady-state  compressible  flow  of  gas  in  pipes. 


9b 


jb  ,no  ^  u  f  o  2  6  nf-fijdo  oJ  '.sbTo  n 


. si rn  rl 

JicrJs  1  "to  b  r  6  arf  ^  ri?  f  (Sf:-A)  bne  (C£-A)  2fforJ6uo3 

26  n93JHw9*T  9d  ysm  (08-A)  HguoifU  (?S-A) 


{’a.U.^.O.^.^.^.J  q.„T)„0^  =  i-b- 


bne 


(K-A) 


,T.  9,  T)«t>t  « 


-  9fiub  eqrq  S(\i  nsrlw  taf  is  At  «92bd  won  oilrosqe  ,  to3 
2sg  bne  (Qp  eJfit  won  t(fctO,*J)  eosq2  ni  ftofJr2oq  bns  2norz 
Ao-A)  brie  (£t-A)  2ftor:f6Up3  tnsvrQ  stb  (  rtff  9,  T)  29 mac o-tq 

c  SOUbS't 

qb 


APPENDIX  B 


DEVELOPMENT  OF  MATHEMATICAL  EXPRESSIONS  FOR  THE  PARTIAL 
DERIVATIVES  OF  Z  WITH  RESPECT  TO  Tp  AND  Pf 


The  reduced  Benedi ct-Webb-Rubin  equation  of  state  has 
the  form 


Z 


A  n  A  o  A  r  O  Py, 

+  (A1  +  T7  +  +  (A4  +  ?7Pr  +  A5A6 

p2 

+  A  7  — r-j  (1  +  AgP^)  exp  (  -AgP^)  (B-1) 

Tr 


where 


Z 


RT  pRT 


(B-2) 


At  the  critical  point 


P  M 

zc  ■  p^T-  (B-3> 

c  c 

Dividing  Equation  (B-2)  by  Equation  (B-3)  and  rearranging 
give 
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Z  =  Z 


c  p  T 
Mr  r 


( B  -  4 ) 


Combining  Equations  (B-l)  and  (B-4)  produces  another 
form  of  reduced  Bendi ct-Webb-Rubi n  equation  of  state 


r  Z 


rTr  +  <AlTr  +  A2  +  +  <Vr  +  Vp? 

r 


+  A5A6pr  +  A7  "^T  ^  +  A8pr^ea?p^"A8pr^ 

Tr 


(  B  -  5  ) 


az 


3  Z 

Equation  (B-4)  was  employed  to  derive  (tt- )p  anc* 

r  r 


(|p— )T  by  partial  differentiation  as  follows: 
Kr  r 


Z  P 

/  9  Z  \  _  9  /  c  rx 

3Tr  Pr  '  3Tr  ‘VrV 


ZcPr 


1  1 


3p, 

(?T-) 


o  T  2  T  d  2  •  3VPr 
Hr  r  r  r 


vv 

pr^r 


.  J_  .  J_  (ffll) 

T  p  V3T  ;P 
r  Kr  r  r 


(|f->p  ■  1 

r  r 


1  1  3pr 

T  pT  ^TT~^P 

r  Kr  r  r 


( B-6 ) 


% 
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And 


(AL.) 

3Pr  Tr 


3 P  VT  p  '  T 
r  rKr  r 


3pr 

3P 


ZcPr 

p77 


1  ,3pr 
-  Pr  (3Pr 


) 


l3F^T 
r  r 


1 '  3pr 

—  ( — -1 
D  ^  3P  'T 
r  r  r 


( B-7 ) 


Equations  (B-6)  and  (B-7)  imply  that  the  partial  deriva¬ 
tives  of  Z  with  respect  to  Tr  and  can  be  evaluated  once 
the  partial  derivatives  of  pr  with  respect  to  Tr  and  Pr  are 
known . 


3pv- 

Now  ( itd  ) 


1 


ap. 


3Pr'Tr  Z7T 
8pr  Tr 


;  provided  (-g^-)T  f  0  (B-8) 


r  r 


3p 

( _ -) 

'p 


r 

r  ’  r 


3P  3P 


( B  -  9 ) 


67  b  Tfir^-sq  siU  \  qm  (X-8)  bn;  (8  i- /  2nof^6Up3 

T  o i  ^osqae-t  rIJrw  X  ’to  asvr i 


josqaat  rlJh  q  zsvtJfiv  .  bb  rerJveq  sriJ 


*i  "i  A  f  *•( 1 


„r  „  V 
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3Pr  3P 

(gTf — ) 0  and  (~—  )t  were  obtained  by  performing  partial 
r  r  pr  r 

differentiation  on  Equation  (B-5) 


3P 

( FT^p 

r  r 


pr  +  ^A1 


2A 


T. 


3x  2  ,  .  3 

-j)pr  + 


pr  /  *i  i  n  2 


"  2A7  — 3  (1  +  Agpr)exp(-Agpp 
Tr 


(B-10) 


9Pr 
r  r 


Tr  +  2<AlTr  +  A2  +  7tK 


+  3<A4Tr  +  A5>Pr  +  6A5A6pr 


+  A 7  -^2*  ea;p(-AgP^)  ( 3  +  3AgP^ 
Tr 


2Agp^) 


(B-ll) 


Substituting  Equation  (B-ll)  into  Equation  (B-8)  gives 


no 


9pr 

^3P~^T 
r  r 


Tr  +  2<AlTr  +  A2  +  A>pr 

r 


+  3<A4Tr  +  A5>pr  +  6A5A6pr 


+  —^2  exp(-A gp2)(3  +  3Agp2  -  2AgpA) 


(B-12) 


Substituting  Equations  (B-10)  and  (B-ll)  into  Equation 
(B-9)  gives 


3pr 

(yH  p 

r  r 


2  A«  «  o 

pr  +  (A1  -  7T)pr  +  A4pr 


-  2A7  — j  (1  +  Agpr)exp(-Agpr) 

- 

A. 


Tr  +  2<AlTr  +  A2  +  7T)p: 


+  3(Vr  +  A5)pr  +  6A5A6pr 


+  A 


7  T  2 
r 


(B-13) 


ea?p(-Agp2)(3  +  3Agp2  -  2Agp A) 


irJti*  M  zdu2 


•'  W*  +  '«<  \  ^  ,1*A)  + 


(!%4aS  -  >gA£  +  e)(^q8A-)q^  ^  + 


Ill 


Equations  ( B - 6 )  and  (B-7)  with  the  aid  of  Equations 

3  7 

( B - 1 2 )  and  ( B - 1  3 )  made  it  possible  to  estimate  (tt- )p  and 

3Z  r  r 

(jp — )T  at  T  and  P  analytically. 


'  < 


APPENDIX  C 


DEVELOPMENT  OF  A  MATHEMATICAL  EXPRESSION  FOR 
HEAT  CAPACITY  DEPARTURE 


The  reduced  form  of  the  Benedi c t-Webb- Rub i n  equation 
of  state  is 


P 


r 


prTr  +  <AlTr  +  A2  +  JWr  +  <Vr  +  Vp3r 

r 

p3 

+  A5A6pr  +  A7  "^T  ^  +  Aspr^ea;p^"A8pr^ 

^r 


(C-l) 


The  heat  capacity  departure  from  the  ideal  state  is 

def i ned  as  ( Ctt  -  C° )  . 

P  P 

For  a  gas  in  the  ideal  state 

c°  -  c°  +  R  (  C- 2 ) 

Thus 

(Cp  -  cp>  =  cp  -  (cv  +  R)  (C‘3) 

By  adding  and  subtracting  Cy  from  the  left  hand  side  of 
Equation  (C-3)  and  rearranging,  one  may  get 
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(S  -  9  =  (CP  -  cv> +  (cv  -  9  - R 


( C-4) 


where  C  =  isobaric  heat  caoacity  at  pressure  P  and  tempera- 


r  ture  T,  Btu/lb  mole.°R. 

o 

C  =  isobaric  heat  capacity  at  zero  pressure  and 
p  temperature  T,  Btu/lb  mole.°R. 

Cv  =  cons  tan t-vol ume  heat  capacity  at  pressure  P  and 
temperature  T,  Btu/lb  mole.°R. 

C°  =  cons tant-vol ume  specific  heat  capacity  at  zero 
pressure  or  at  infinite  volume  and  temperature 
T ,  B t u / 1 b  mol e . 0 R. 

R  =  universal  gas  constant,  1.98654  Btu/lb  mole.°R. 

Expressions  for  (C^  -  Cv)  and  (Cv  -  C°)  may  be  obtained  from 
Equation  (C-l)  by  applying  the  following  thermodynamic 
relationships  ( 176) : 


_ V) 

9  V  '  T 


3Cv 


(C-5) 


<Cp  -  Cv>  = 


(C-6) 


Equations  (C-5)  and  (C-6)  may  be  rewritten  as: 


(C-7) 


( C-8) 


err 


-B'tsqnteJ  bne  8  9iU229iq  *6  vd/06q60  deed  oM&dozr  -  .0 

.K°.9rorrf  df\ui  l  ,T  e-u/d 

bne  9iu229*fq  ot^s  d6  ^j-fDBqso  d*>sd  orTedoar  =  *0 

fi.sfom  df\ud8  ,T  9tuj  t  •  •.,  j,-:9d 

b  is  9  9*i  u  v.  2  9't  q  26  '{drofiqsD  dfced  amu  Tov -3  ns  3  ?,nci  -  3 

.^°.9ro:iT  dr\ud8  tT  9TU dfi'tqqnisd 

CHS!  d6  y,dfD6q6D  J69fl  0  3  I"  D  9q  8  9fnU  TOV  -d  n  6d  2  ftOD  =  °3 
squJfiq9qin9i  brse  9jnurov  9drnrdnr  d6  to  a'ojaa  'q 

.H° . 9  Torn  dr \ud8  t  T 

.fl°.9[om  df\ud£  £<?d8G.f  ,  dnedanoo  a&g  ffia'isvrnu  = 


mont  bsn  red  do  9d  \sm  (°3  -  r0)  bnfi  (y3  -  3)  qod  anoraaeqqxH 

o  rmsnYborrrorid  gniworrod  odd  gniy,rqq6  (r*3)  rord&u^3 

:(dKT)  zq r rtano rdsf sq 


(6-3) 


0-3) 


S’ 6  w  S,S6 


(V°  -  q3) 


:as  nsJJjiwsT  sd  0-0)  bus  O-D!  sriofdf.i;p3 


(V-D) 


(8-0) 
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By  performing  partial  differentiation  on  Equation 
(C-l),  the  following  expressions  were  obtained: 


3P 

(TT>p 

r  Kr 


Pr  +  (A 


2A 


1  -TT'^r 


3x  2  ,  .  3 

)P.  +  A4Pr 


-  2A 7  -^3  (1  +  Agp2 ) exp ( " AgP^ ) 
^r 


( C-9 ) 


32P 
( - 'o) 


1 


3Tr2  pr  Zc 


6A3  rr  +  6A7  rV  (1  +  Vr^^'V^ 
Tr  Tr 


( C-l 0) 


3P  i 

r  r  c 


Tr  +  2(A1Tr  +  A2  +  -^)pr 


+  3<A4Tr  +  A5)pr  +  6A5A6pr 


+  A 7  -^2  e^p(-Agp2)(3  +  3Agp2 
Tr 


2A8pr ) 


(C-l 1 ) 


Substituting  Equations  (C-9),  and  (C-l  1 )  into  Equation  (C-7) 
gives 


)  I  -t  £  u  •  j  ro  no  r  Je  rb  fe-jieq  r-n  rnno'trt9q  y3 

:b9n;6}do  s^isw  anor  ?a-3iqx9  gniwof  o'*  bd)  t(f-3) 


(6-3) 


+  V  <f 


A)  +  *!  Q 


'f-/  "  "  '  '  Qr;;  +  )  ■;  v  ■ 


(or-3) 


{^Q8A-)qx9(^q8A  +  f)  j2-  ^A3  +  j'-  fAc 


) 


c-  -  r  6 


(rr-o) 


($c,.  AS 


+  sA  +  ,TrA>s  +  -,T 


;;qdAaA3  +  >(aA  +  ,T,A)e  ♦ 


^QgAC  +  C) (^dgA-)qxs  j 

•t 1 


r  r  -  3 )  brte  «.  ( Q  -  3 )  a  n  o  i  ;t  e  ij  p  3  gnMudriaduE 

29V  rg 
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-  2 


+  (A- 


2A 


|)Pr  +  A4P2r 


-  2A7  -^3  (1  +  Agp2 ) exp ( ~ AgP^ ) 
Tr 


Tr  +  2<AlTr  +  A2  + 

r 

+  3(A4Tr  +  Vpr  +  6A5A6pr 
p2 

+  A j  — C~2  exp(-AQp^)  (3  +  3Agp2  -  2AgpJ) 
Tr 


(C-12) 


Substituting  Equation  ( C - 1 0 )  into  Equation  ( C - 7 )  gives 


ac 

(_ V) 
l3p  ' T 

pr  r 


=  -  R 


6A- 


T. 


3  +  6A7 


T. 


-j  0  +  Agpr)exp(-Agpr) 


(C-13) 


Boundary  Conditions: 

At 

By  separating  variables  in  Equation  (C-13),  integrating,  and 
applying  the  above  Boundary  Conditions,  one  may  obtain: 
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// (8Cv>T  =  -  R 


,P-  6A7  p  p  ~ 

/  ( — 3)7  dp  +  /  (6Ay  — 3  exp(- ArP  ))T  dp 

0  T  J  r  r  0  '  T  J  y  r  V  1 

r  r 


P  P 

r/^„  „  Hr 


+  /  (6A7A8  fT  eaP('A8pr)  )T| 


dp 


(Cy  -  Cv)  =  -  R[X  +  Y  +  Q] 


(C-14) 


where 


X  =  6A 


3  T  3 
r 


(C-15) 


3A?  exp(-AgPr)  3A?  1 

Y - 3 - +  “AT  tt 


A8Tr 


8  T. 


(C-16) 


^  =  "3A7  7T  ^P('A8Pp  -  -fl- 
Tr  8 


2  3A7  ejjp(-AgPr)  3A7 


( C- 1 7 ) 


or 


CCv-Cy)  =  R 


Pr  6A7 

-  6A3  — 3  -  - 3-  +  ( 

0  T  0  AT 

'r  0  r 


6A 


7  +  3A-7  -^-)exp(-Agp^) 


A  t  3  7  j  3 

a8  r  'r 


(  C  -  1  8 ) 


Substituting  Equations  ( C - 1 2 )  and  ( C -  1 8 )  into  Equation 
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(C-4)  results  in  the  final  form  of  heat  capacity  departure 
equation,  which  may  be  written  as: 


(S-9 


=  R 


-  1  -  6A 


pr  6A7 


3  t  ^  'at3 
‘r  o  r 


+  (- 


6A 


Vr 


^3  +  3A7  ^-j)e«p(-AgP2) 


T 


-i  2 


2A. 


1  +  ^A1  -  I"4)pr  +  A4pr 
r 


-  2Ay  —^3  (1  +  Agp2 ) exp ( " Agp2 ) 
^"r 


A. 


Tr  +  2(AlTr  +  A2  +  ^T)pr 


+  3<A4Tr  +  Vpr  +  6Wr 


+  Ay  -^2  exp(- Agp3)(3  +  3Agp3 
Tr 


2A8pr ) 


(C-19) 


APPENDIX  D 


THE  HYDRAULIC  FRICTION  FACTOR  IN  BOTH  THE  TRANSITION 

AND  THE  ROUGH  PIPE  REGIONS 

Transition  Region 

This  region  is  characteri zed  by  4000  <  Rg  <  ^200 
The  hydraulic  friction  factor  is  given  by  the  Colebrook 
transition  law  (177). 


f 


1o^4t) 


2.51 


(D-l) 


where  f  =  Moody's  friction  factor. 

e  =  absolute  roughness,  inch, 
d  =  I . D .  of  tubing,  inch. 

Re  =  Reynolds  number. 

By  taking  the  square  root  of  both  sides  of  Equation  (D-l) 
and  by  letting  x  =  /  f  ,  one  may  obtain 


x 


2 


1 


log( 


e/6  ,  2 . 5 1 \ 

3.7  R  x  1 
e 


( D-2 ) 


In  order  to  solve  for  x  in  Equation  (D-2),  the  Newton- 
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Raphson  iterative  method  was  used.  The  iterative  equation 
has  the  following  form: 


gUJ 

x.  +  ,  =  x.  -  —  ;  providing  g  (x.)  f  0 

9  (*}) 


(D-3) 


where  g(x.)  is  obtained  from  Equation  (D-2)  as: 


9(xi)  =  1  +  2xi  +  |i|i)  ( D - 4 ) 

e  i 


Thus  , 


g"(xi ) 


2  log{^$-  + 


2.51 

Rexi 


Rexi( 


5.02 

(e/d) 

3.7 


log 


2.51 

Rexi 


(D-5) 


Combining  Equations  (D-3),  (D-4),  and  (D-5)  provides  the 
final  form  of  the  iterative  equation  for  x, 


1  +  2x 


+  2^51, 

Rex1  -1 


2  log( 


(e/d)  ,  2 . 51 x 

3.7  R  x . 

e  i 


5.02  log  e 


Rexi< 


(e/d)  +  2-51 
3.7  Rex. 


) 


(D-6) 


Once  Equation  (D-6)  is  satisfied  after  n  iterations 
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The  Fully  Rough  Pipe  Region 

3200 

Flow  of  fluid  falls  in  this  region  when  R0  >  ^ £yd ^  . 
The  Van-Karman  equation  was  employed  to  estimate  the  hydrau¬ 
lic  friction  factor  in  this  region. 


f  = 


2 


(D-7) 


.  •  <  ft  nsriw  norgat  sNi  nf  affs^  bu  \  woM 

-1/6'  b^ri  9riJ  9^6m’r^29  o*  b9^ofqrns  aew  rtorisups  r<6nn6>l-n6V  9dT 

.  norge**  2 r ni  norj-DM^  off 


( ^  -0 ) 


APPENDIX  E 


AN  ALGORITHM  FOR  SOLVING  THE  SYSTEM  OF  DIFFERENTIAL 
EQUATIONS  IN  MAIN  MODELS  AA  AND  BB 


The  system  of  differential  equations 

that  describes 

the  flow  problem  is  as  follows 


dpr 

^  =  f1(Tr,Pr,R*,£/d) 

(E-l) 

dT 

3T1  =  f2(WR*,£/d) 

(E-2) 

dR*  -  0 
dL  "  0 

(E-3) 

The  Boundary  Conditions  are: 


At  L  =  0  :  T  -  T(l.0)  -  T, 

L  =  L*  :  T  =  =  =  T2  and  P 

P ( L= L * )  '  P2 

The  modified  Runge-Kutta-Merson  algorithm  as  suggested 
by  Niesse  (178),  Chai  and  Burgin  (179)  was  utilized  to  solve 
this  system  of  differential  equations  as  follows: 

(1)  Let  (e/d)  =  0.0006/d  (for  new  tubing) 

R*  =  7000  ;  1000  <  R*  <  7000 
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.  .  .  h  •  n 
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(2) 


For  a  trial  step  size  (h^),  compute  kp  and  ky  ,  i = 1 , 5 


X  VTr  ’ Pr  ’R*’e/d>  i 


n  n 


kT,  =  X  f2<Tr  •Prn’R*’e/d) 

n  n 


kp  =  -f  fi (T  +  kT  ,P  +  kp  ,R*,e/d)  ; 
p2  J  1  rn  'l  rn  P1 


;T  =  -f  f  ,  ( T  +  kT  .P,,  +  kD  , R* , e/d ) 

'2  J  rn  'l  rn  *1 


h  kT  kT  kp  kp 

kp  -  f,(T  +  ~r  +  ~r  .pr  +  ~y-  +  ,R*,£/d) 

r  3  n  n 


kj  =  -y-  f,(T  +  -p-  +  -p-  ,P  +  -j!  +  — ,R*,E/d) 
3  r  n  r  n 


k  =  —  f(T  +  —  k  +  —  k  P  +  —  k 

P4  3  ll \  8  T1  8  V  rn  8  P1 


+  |  kp  ,R* , e/d) 


kI4  ’  3k  f2<Trn  +  8  kT,  +  8  kT3>Prn  +  8  kP, 


+  I  kp  ,R*  ,e/d) 
3 
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kp  =  4  ML  +  |  kT  -  I  kT  +  6kT  ,P 
k5  j  '  rn  ^  1 1  ^  1 3  '4  rn 


+  I  kP,  'I  kP3  +  6kp  ,R*,e/d) 


^  U  "5  Q 

kT  =  ~ T-  fn(T  +  «•  kT  -  -K  kT  +  6 kT  ,P 
T5  ^  2  rn  2  T1  2  T3  T4  rn 


+  |  kD  -  %  kD  +  6kD  , R* , e/d ) 
L  K1  2  k3  p4 


and  calculate 


^rn+l  =  Prn  +  2  Kpl  '  2  kp3  +  6kp4 


\+l  =  T|"n  +  2  “T1  "  2  kp3  +  'S 


(3)  Compute  the  estimated  truncation  error 


EP  ■  pc 


kp  2  + 


EI  "  Tc 


l<  -  —  k  +  4k  -  _ — 

kT]  2  kT3  qkT4  2 


(4)  Let  6 i  =  tolerance  in  predicted  pressure  P. 

62  =  tolerance  in  predicted  temperature  T. 

If  either  Ep  >  6-,  or  ET  >  $2,  halve  the  trial  step  size 
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and  start  oyer  from  step  (2)  until  Ep  1  6^  andEy<  6^. 


(5)  Compute  the  corrected  value  of  P  and  T  : 

rn+l  rn+l 


P  (corrected)  =  P 
rn+l  rn+l 


T  (corrected)  =  T  -  r- 
rn+l  rn+l  c 


This  correction  raises  the  accuracy  to  the  fourth  order. 


6 1 

(6)  If  Ep  <  jg-  and  Ey  <  yg-  ,  double  the  trial  step  size  for 
the  next  integration  step.  If  not,  leave  the  trial  step 
size  as  is  for  the  next  integration  interval,  h^  . 


k 

If 

L*  - 

l 

i  =  1 

hi 

=  0  , 

go  to 

step  (8) 

— 

k 

— 

If 

L*  - 

I 

hi 

i  hk+l 

,  go 

to  step 

(2). 

— 

i  =1 

— 

k 

k 

If 

L*  - 

I 

i  =1 

hi 

A 

zr 

+ 

,  1  et 

hk+l  = 

L*  - 

I  h. 

i  =  1  J 

and  go  to  step  ( 2 ) . 


(8)  At  this  stage  T^l=0j  and  P(l  =  0)  are  estimated. 

1 f  I(L=o)  =  T1  ’  then  P1  =  P(L=0^  *  ST0P* 


' 


If  T(l  =  0)  +  T-j  ,  modify  the  value  of  R* 


(9)  Let  k  =  1  and  start  over  from  step  (2) 


For  main  model  BB,  R**  replaces  R*. 


APPENDIX  F 


SOURCE  DATA 


WELL  A 


FIELD :  Dunvegan  WELL  NAME:  Anderson  et  al . 

Dunvegan  6-29 


POOL:  Debolt  LOCATION:  6-29-80-3  W6M 


PERFORATIONS 

:  4748' 

-4785' 

K.B. 

DISTANCE 

:  TO  M. 

P.P. 

:  4753'  K.B 

TUBING  SIZE: 

2.441" 

I.D. 

RESERVOIR  TEMPERATURE:  115°F 

GAS  PROPERTIES: 

G  =  0.6402 

pPc 

669. 

9  PSIA 

PTC  ■ 

367. 

2  °R 

mol.%  co2  = 

0.56 

MOL.% 

N2 

=  1.05 

MOL.% 

H2S 

=  0.00 

MOL . %  C1 

90.33 

MOL.% 

C2 

=  4.04 

MOL.% 

C3 

=  1  .95 

MOL . %  i -C 4  = 

0.31 

MOL.% 

n-C4 

=  0.70 

MOL.% 

^5 

=  0.20 

MOL . %  n-C5  = 

0.25 

MOL.% 

r  »  S 
L6 

=  0.22 

MOL.% 

r  .s 

=  0.25 

MOL.%  C+ 

0.07 

MOL.% 

He 

=  0.07 

Density  of  Condensate  =  0.777  gm/c.c. 


WELL  OPERATOR:  Anderson  Exploration  Limited 
TYPE  OF  TEST:  Modified  Isochronal 
DATE  OF  TEST:  April  5,  1971 
GAS  PRODUCED  THROUGH:  Tubing 
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WELL  A  (CONTINUED) 
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WELL  B 

FIELD:  Bigstone  WELL  NAME:  Pan  Am.  HB,G-2 

POOL :  D-3A  LOCATION:  2-25-60-22  W5M 

PERFORATIONS:  n006,-11052'  K.B.  DISTANCE  TO  M.P.P. :  11029'  K.B. 
TUBING  SIZE:  2.992"  I.D.  RESERVOIR  TEMPERATURE:  246°F 

GAS  PROPERTIES: 

G  =  0.6997  P  =  801.2  PSIA  T  =  410.9  °R 

pc  pc 

MOL.%  C02  =  2.88  MOL . %  N2  =  4.56  MOL.%  H2S  =  18.89 
(Complete  gas  analysis  at  the  time  of  test  is  not  available.) 

WELL  OPERATOR:  Amoco  Canada  Petroleum  Company  Limited 

TYPE  OF  TEST:  Multi-Point 

DATE  OF  TEST:  September  20,  1970 

GAS  PRODUCED  THROUGH:  Tubing 


WELL  B  (CONTINUED) 
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WELL  C 


FIELD:  Bigsto 

POOL:  D-3A 

PERFORATIONS: 

TUBING  SIZE:  2 

GAS  PROPERTIES 
G  =  0.702 

mi.%  co2  =  2. 

(Complete  gas 
WELL  OPERATOR: 

TYPE  OF  TEST: 

DATE  OF  TEST: 


ne  WELL  NAME:  Pan  Am.  HB.C-l 

LOCATION:  13-35-60-22  W5M 

1 0957 '-1  0974'  K.B.  DISTANCE  TO  M.P.P.  :  10965'  K .  B  . 
.992"  I.D.  RESERVOIR  TEMPERATURE:  243°F 

P  =  798.3  PSIA  T  =  409.5  °R 
pc  pc 

78  MOL . %  N2  =  4.56  MOL . %  H2S  =  18.50 
analysis  at  the  time  of  test  is  not  available.) 

Amoco  Canada  Petroleum  Company  Limited 

Single-Point 

October  10,  1972 


GAS  PRODUCED  THROUGH:  Tubing 


. 
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WELL  D 

FIELD:  Bigstone  WELL  NAME:  Pan  Am.  HB,G-2 

POOL:  D-3A  LOCATION:  2-25-60-22  W5M 

PERFORATIONS:  1  1006,-n052'  K.B.  DISTANCE  TO  M.P.P.  :  1  1029'  K .  B . 

TUBING  SIZE:  2.992"  I.D.  RESERVOIR  TEMPERATURE:  243°F 

GAS  PROPERTIES: 

G  =  0.702  P  =  798.3  PSIA  T  =  409.5  °R 

pc  pc 

MOL . %  C02  =  2.78  MOL . %  N2  =  3.52  MOL.%  H2S  =  18.50 
(Complete  gas  analysis  at  the  time  of  test  is  not  available.) 

WELL  OPERATOR;  Amoco  Canada  Petroleum  Company  Limited 

TYPE  OF  TEST:  Single-Point 

DATE  OF  TEST:  October  10,  1972 

GAS  PRODUCED  THROUGH:  Tubing 
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APPENDIX  G 


COPIES  OF  THE  FORTRAN  COMPUTER  PROGRAMS 
FOR  MAIN  MODELS  AA  AND  BB 

C  *************************************************** 
C ************ *******  MAIN  MODEL  AA  **************** 

C*************************************************** 


C **  THIS  MAIN  MODEL, 'AA',  REQUIRES  THE  FOLLOWING 

C**  FUNCTION  AND  SUBROUTINE  PROGRAMS  : 

C**  F, SEARCH, AMEIN,KASEM,BWR, LEE, MOODY, 

C**  QUON  &  MISIC  .  THESE  PROGRAMS  ARE  LISTED 

C**  IN  PAGES  TO  FOLLOW  . 

C**  HOW  TO  USE  THIS  COMPUTER  PROGRAM 


C 

C 

C 

C 

C 

C 

c 

c 


ITYPE  IS  A  FLAG  INTEGER 

ITYPE=0  :  INJECTION  PROBLEM 
ITYPE=1  :  PRODUCTION  PROBLEM 
M  IS  A  FLAG  INTEGER 

M=1  :  SUPPLY  GAS  ANALYSIS 

MMM  IS  A  FLAG  INTEGER 

MMM=0  :  ONE-PHASE  FLOW 
MMM= 1  :  TWO-PHASE  FLOW 


C 

c 

c 

c 

c 

c 

c 


COMBINATIONS  OF  THE  FLAG  INTEGERS  ARE  AS  FOLLOWS: 
ITYPE  M  MMM  DESCRIPTION 

1  1  1  TWO-PHASE  FLOW  PRODUCTION 

PROBLEM.  SUPPLY  GAS  COMPOSITION 
AND  SG  OF  FLOWING  STREAM 
1  1  0  ONE-PHASE  FLOW  PRODUCTION 

PROBLEM.  SUPPLY  GAS  COMPOSITION 


C  GASUC  IS  GAS  UNIVERSAL  CONSTANT 

C  GC  IS  CONVERSION  FACTOR 

C  AJ  IS  CONVERSION  FACTOR 

C  SLOPE  IS  SLOPE  OF  PIPE  SEGMENTS 


c**  THE  FOLLOWING  TERMS  BELONG  TO  PURE  GASES 
C  YMOLF  IS  MOLE  FRACTION 

C  YMW  IS  MOLECULAR  WEIGHT 

C  YMW1  IS  MOLECULAR  WEIGHT 

C  T EMC  IS  CRITICAL  TEMPERATURE  IN  DEG.  RANKIN 

C  TEMC1  IS  CRITICAL  TEMPERATURE  IN  DEG.  RANKIN 

C  PREC  IS  CRITICAL  PRESSURE  IN  PSIA 

C  PREC1  IS  CRITICAL  PRESSURE  IN  PSIA 

C  VOLC  IS  CRITICAL  VOLUME  IN  CU.  FT/LB 

C  VOLC1  IS  CRITICAL  VOLUME  IN  CU.  FT/LB 

C  TNEOI  IS  NORMAL  BOILING  POINT  IN  DEG.  R 

C  TNBOI  IS  NORMAL  BOILING  POINT  IN  DEG.  R 

C  SEGM1  IS  COLLISION  DIAMETER  IN  ANGESTROMS 

C  SEGMA  IS  COLLISION  DIAMETER  IN  ANGESTROMS 


-p-  UJ  tv)  -*  t  UJ  M  ->  -pr  (_U  NJ 
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C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


EOK  IS  A  PARAMETER  IN  LENNARD- JONES  POTENTIAL, 
IT  IS  (EO  /K  )  IN  THE  TEXT 
E OK 1  IS  THE  SAME  AS  EOK 

CPOC  IS  THE  FOUR  COEFFICIENTS  OF  THE  IDEAL 
HEAT  CAPACITY  EQ.  OF  PURE  GASES 
(KOBE  ET  AL.  ) 

CPOC1  IS  THE  SAME  AS  CPOC 


SG 

IS 

PTEMC 

IS 

PPREC 

IS 

DIAM 

IS 

DEPTH 

IS 

TSOD 

IS 

TSD 

IS 

TDD 

IS 

PSD 

IS 

Q 

IS 

GAS  GRAVITY 

PSEUDOCRITICAL  TEMP.  IN  DEG.  RANKIN 
PSEUDOCRITICAL  PRESS.  IN  PSIA 
TUBING  DIAMETER  IN  INCHES 
LENGTH  OF  TUBING  IN  FEET 
GROUND  SURFACE  TEMP.  IN  DEG.  F 
TUBING  HEAD  TEMP.  IN  DEG.  F 
BOTTOM-HOLE  TEMP.  IN  DEG.  F 
TUBING  HEAD  PRESSURE  IN  PSIA 
TOTAL  FLOW  RATE  IN  MSCF/DA Y 


IMPLICIT  REAL*8  (A-H,0-Z) 

EXTERNAL  F 

DIMENSION  VOLC1  (15)  ,TNB01  (15)  ,SEGM1  (15)  ,E0K1  (15) 
DIMENSION  YMW  1  (15)  ,TEMC1  (15)  ,PREC1  (15) 

DIMENSION  CPOC1  (15,4) 

COMMON  /BL1/GASUC,GC,AJ, SLOPE/BL 2/1 TYPE 
COMMON  /BL1C/TNBOI (15) /BL4/SEGMA  (15)  ,EOK(15) 
COMMON  /BL5/YMW  (15)  ,TC  (1  5)  ,  ZC  (1  5)  , GAMA  (15) 

COMMON  /BL8/CPOC  (15,4) 

COMMON  /BL11/TEMC  (15)  ,PREC  (15)  ,VOLC  (15) 

COMMON  /BL12/PTEMC,ZT EMC, PPREC, Z PREC/BL 18/M 
COMMON  /BL13/DIAM/BL1 4/SG/BL9/YMOLF (15) 

COMMON  /BL1 5/Q,TSD, PSD, TDD, DEPT H,TSOD,DL 
COMMON  /BL30/VMASS 
COMMON  /BL50/MMM 

DATA  YMW 1/C .  440 1 3D  02,C.  28013D  02,0.34C76D  02, 

1  0.16043D  02,0. 3G07GD  02,0.44C97D  02, 

0  2,0 . 581 24D  02,0.72151D  02, 
02,0.86178D  02, 0. 66 1 78D  02, 
03,0. 10021D  03, 0. 1 1 423D  03/ 
03,0. 22727D  03,0.67237D  03, 
03,0 . 54976D  03,0.66568D  03, 
03,0.76532D  03,0.82877D  03, 
03,0. 89 550D  03,0.91337D  03, 


DATA 


0. 58124D 
0.72151D 
C. 10021 D 
TEMC1/0. 54757D 
0. 34304D 
0.73465D 
0. 84537D 
0. 95467D 


DATA  PR EC  1/0 . 107 ID  04,0.493CD  03,G.13C6D  04, 
0.6678D  03,0. 7078D  03,0.6163D  03, 
0.5291D  03,0. 55C7D  03,0.4904D  03, 
0.4886D  03,0. 4366D  03,0.4369D  03, 
C.  3965D  03 ,0 . 3968D  03,C.3606D  03/ 
DATA  VOLC1/0. 342D-01 ,0.514D-01 ,0 . 459D-01 , 

1  0. 991D-01,0.788D-01,0.737D-01, 

2  0.724 D- 01, C.7C2D-C1,Q.679D-01, 


• 

■iz 

B  ■  ) 

: 

- 

. 

4 


C  U)  to  -»  -P  OJ  K>  -C  UJ  to 


136 


3  0.675D-01, 

4  0 • 673D-0 1 , 
DATA  TNB01/0. 35037D  0 

1  G.2C098D 

C. 47057D 
0. 55659D 
0. 65376D 

DATA  SEGM1/0. 3996D 
C. 3758D 
C. 5278D 
0. 5784D 
0. 7C15D 
DATA  EOK1/0.19C0D 
0. 1486D 
C.33G1D 
0 . 34 1 1D 
C. 2742D 
DATA  CPOC1/0. 5316D 
0. 4750D 
-0. 1890D 
C. 1618D 
-C. 26069D 


A 

B 

C 

D 

E 

F 

G 

H 

I 

J 

K 

L 

M 

0 

P 

Q 

R 

S 

T 


0 

0 

c 

0 

01 

01 

01 

01 

01 

03 

03 

03 

03 

03 

01 

01 

01 

01 

0 


0. 793611  ID-02, -0 
0. 6666666D-02,  0 
0. 5519999D-01,  0 
0.6027777D-01,  0 
0. 961 1055D-Q1 ,  0 
-C. 2580864D-05,  0 
0. 93518510-06,-0 
-0. 1695987D-04,-Q 
-0. 1655864D-04,-0 
-0.  30  89506D-04 , -C 
0. 3058984D-09,-0 
-0.4509602D-09,  0 
0 • 20  4389  5D-G8 ,  0 
G. 1731824D-08,  0 
C. 3834362D-08,  0 
GC=0 . 32 1 74D  02 
GASUC=C. 1 545D  04 
AJ=0.77816D  03 
SLOPE=0 • ID  01 


0*6810-01,0, 68 
0. 691D-0 1 ,0.69 
3,0.1 3927D  03, 
3,0. 3321 9D  03, 
3,0. 49077D  03, 
3, 0 . 600 1 4D  03, 
3,0. 66 884D  03, 
, 0 . 379 8D  01 ,C. 
, 0 . 4443D  01,0. 
,0. 4687D  01,0. 
,0. 6256D  01  ,0. 
,0.7C04D  01 ,0. 
,0 . 71 4CD  02,0. 
,0.21 57D  03,0. 
,0. 531 4D  03,0. 
,0 . 2472D  03,0. 
,0.2882D  03,0. 
,  0.  6903D  01  , 

,  0.1648D  01,- 
,  0.9450D  00,- 
,-0. 9081D  00, 
1,  0.31899D  01 
.20  8  5GC0D-0  3 , 

•  2291 1 1  1D-01 , 

• 4929444D-0 1 , 

• 7843494D-01 , 
.831 56C5D-0 1 , 
.5956790D-06, 

•  4722222D-0  5,- 
. 1351 851D-Q4,- 

•  23969 10D-0  4 ,- 
. 23731 48D-04 ,- 
. 1 176440D-09,- 
• 2983539D-09, 

.  14334  70D-0  8 , 

•  2  84  8422D-0  8 , 
.261 4881 D-0  8, 


8D-01 , 

OD-0 1  / 

0.38307D  03, 
0.41600D  03, 
0.54179D  03, 
0.61539D  03, 
0.71789D  03/ 
3623D  01, 

5118D  01, 

605 1 D  01, 

5949D  01, 

7451 D  01/ 

30 1 1 D  C 3, 

2371 D  03, 

2694D  03, 

3993D  03, 

3200 D  03/ 
Q.7070D  01, 
0.9660D  0C, 
0.2273D  01, 
0.1657D  01, 

,  C.35844D  01, 
0.1 737777D-02, 
0.4043888D-01, 
0 . 6907777D-C  1 , 
0. 7327777D-C  1 , 
0. 9502444D-0  1 , 
0.42C9876D-06, 
0. 1 1 58950D-C 4, 
0.21 90432D-C  4 , 
0. 2112345 D- 04, 
0.  2742592D-0  4 , 
0. 1348936D-C9, 
0. 1 299725D-0  8 , 
0. 2719478D-08, 
0. 2362825D-C  8, 
0. 30 6001 2D- 0  8/ 


DO  5  1=1 , 15 
Y MW  (I)  =  YMW1  (I) 

T EMC  (I)  =TEMC  1  (I) 

PREC  (I)  =P EEC  1  (I) 

VOLC  (I)  =V CLC 1  (I) 
TNBCI (I) =  TNB0 1  (I) 
SEGMA (I) =  SEGM 1  (I) 

EOK (I) =  EO  K 1  (I) 

DO  5  J=1, 4 

5  CPCC  (I,  J)  =CPOCl  (I,  J) 
WRITE  (6,208) 


. 

: 


o  o 
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WRITE  (6,202)  (YMW  (I)  ,  TEMC  (I)  ,PREC  (I)  ,  VOLC  (I)  , 
*TNBOI  (I)  ,  SEGMA  (I)  ,  EOK  (I)  ,1=1  ,  15) 

WRITE  (6,213) 

WRITE  (6,2  05)  (  (CPOC  (I,  J)  ,  J  =  1 , 4)  ,  I  =  1 ,  1  5) 

WRITE  (6,2^4) 

C  READ  3  CARDS.  GAS  COMPOS ITION. EACH  CARD  HAS 

C  THE  MOLE  ERACTIONS  OF  5  COMPONENTS.  THE 

C  ORDER  OF  COMPONENTS  IS  AS  FOLLOWS  : 

C  C02, N2,H2S,C1,C2,C3, I-C4, N-C4,I-C5, N-C5, 

C  I-C6, N-C6 ,I-C7,N-C7, N-C8  . 

C  IF  GAS  COMPOSITION  IS  NOT  AVAILABLE, 

C  SUPPLY  3  BLANK  CARDS  . 

READ  (5, 20  6)  ( YMOLF  (I)  ,1=1,15) 

WRITE  (6,215) 

WRITE(6,2C7)  (YMOLF  (I)  ,1  =  1 ,15) 

READ  1  CARD  (  COMBINATION  OF  FLAG  INTEGERS 
AND  SG  ) 

READ (5, 209)  ITYPE , M , MMM , SG 

C  READ  2  CARDS  (  WELL  PRODUCTION  DATA  ) 

READ (5 , 20  6)  DIAM, DEPTH, TSOD 
READ  (5, 20  6)  Q , TSD, PSD  ,TDD 
IF  (ITYPE)  1,1,2 

1  WRITE  (6,218) 

GO  TO  3 

2  WRITE  (6,219) 

WRITE  (6,208) 

3  WRITE (6,2 16)  DIAM, DEPTH 

WRITE  (6,217)  Q, TSOD, TSD, TDD, PSD 
CALL  AMEIN 

WRITE  (6 ,2  27)  SG , PTEMC , PPREC 
WRITE  (6,225) 

WRITE  (6,226) 

C  CALCULATING  MOLAL  FLOW  RATE 

TR=  (0. 6D  02+0 •  45967D  03) /ZTEMC 
PR=0 . 1 465 D  02/ZPREC 

CALL  BWR  (TR, PR , ITER , ZST, DR , CPDM , DZDTR , DZDPR) 
VMASS=0.55744134D-02*Q/(ZST*DIAM**2) 

C  READ  1  CARD  (  RANGE  OF  R*  ) 

C  XMI N= 1 COO . C 

C  XOUT= 1010.0 

C  XMAX=7CG0 • 0 

C  IF  CONVERGENCE  IS  NOT  ACHIEVED, THE  RANGE 

C  MAY  BE  PUT  BETWEEN  200. C  &  1000.0  OR 

C  7000. C  &  2C0000.0 


" 
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C**  REMARK  **  FOR  SAVING  COMPUTER  TIME  IN  CASE  5* 

C  IS  OUTSIDE  THE  SPECIFIED  RANGE, 

C  ALWAYS  PUT  XOUT=XMIN+1 0. C 

READ  (5,201)  XMIN,XOUT,XMAX 
TEMESP=0. ID  00 

C  MINIMIZING  FUNCTION  F 

CALL  SEARCH  ( F , TEMESP ,X , XMI N , XMA X , XOUT) 

201  FORMAT  (7F10. 4) 

2C2  FORMAT  (10X,7F16. 7) 

2C4  FGRMAT  (4D15. 7) 

2C5  FORMAT  (10X,4D20. 7) 

206  FORMAT  (5F15. 8) 

207  FORMAT  (18X,F20. 5) 

2G8  FORMAT  ('  1  ') 

209  FORMAT (3l5,5X  ,F10. 4) 

211  FORMAT  (10X,5F15. 7) 

213  FORMAT (///////, 29X , • A  1 , 1 9X , ' B • , 1 9X, • C ' , 1 9X, • D « , /) 

214  FORMAT  (//,47X  , 'CPO  (I) =A+B*T+C*T**2+D*T**3  • ) 

215  FORMAT ('1 ',////////, 27X, 'MOLE  FRACTION',/) 

216  FORMAT  (//////, 10X,' TUBING  I.  D. • , 27X , '  = » , FI 0 . 3 , 

A  •  INCHES',/, 

B10X, 'THE  LENGTH  OF  FLOW  STRING  =', 

C  F1C.3,'  FEET') 

217  FORMAT (10X, ' FLOW  RATE ' , 3GX, ' =• ,F 10.  3, 

A  '  MSCF/DAY ' ,/, 

B10X, 'THE  AVE.  TEMP.  OF  THE  GROUND  SURFACE  =', 

C  F1C.3, '  F',/, 

D10X, 'THE  TEMP.  OF  FLOWING  STREAM  AT  SURFACE  =', 

E  F10.3, '  F' ,/, 

F 1 0X , *  THE  TEMP.  OF  THE  PRODUCING  FORMATION  =', 

G  F10.3, '  F« ,/, 

HI  OX , ' THE  WELL  HEAD  TUBING  PRESSURE  =', 

I  F1C.3,'  PSIA',/) 

218  FORMAT (10X, 'INJECTION  PROBLEM') 

219  FCPMAT { 1 0  X , 'PRODUCTION  PROBLEM') 

225  FORMAT (1CX, ' THE  CALCULATED  BOTTOM-HOLE  PRESSURE', 

1  '  IS  THAT  WHEN  THE  ABSOLUTE  DIFFERENCE  BETWEEN', 

2  /,10X, 'EARTH  AND  CALCULATED  FLUID  TEMPERATURES', 

3  '  AT  THE  BOTTOM  OF  WELL  IS  WITHIN  0.1  DEG.  F  .') 

226  FORMAT (/, 10X, 'FIRST  COLUMN  IS  DEPTH  IN  FEET  .',/, 

1  1CX, 'SECOND  COL.  IS  EARTH  TEMP.  IN  DEG.  F  .',/, 

2  10X, 'THIRD  COL.  IS  CALC.  FLUID  TEMP.  IN  DEG.  F  .' 

3  ,/, 10X, ' FOURTH  COL.  IS  CALC.  PRESS.  IN  PSIA  .', 

4  /,10X, 'FIFTH  COL.  IS  R*  . ' ,/) 

227  FORMAT  (10X, 'GAS  GRAVITY,  PSEUDO  REDUCED  TEMP.,', 

1  'AND  PSEUDO  REDUCED  PRESS.  ARE  : ' , /, 1 0X, 3F1 0. 4 ) 

STOP 

END 


' 

■ 

' 


139 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


FUNCTION  F  (X) 


THIS  FUNCTION  SUBROUTINE  IS  TO  BE 
USED  WITH  MAIN  MODEL  AA  ONLY 


F 

IS 

E 

IS 

DIAM 

IS 

DL 

IS 

TGEOG 

IS 

PEE 

IS 

TEM 

IS 

DENS 

IS 

REN 

IS 

PRAND 

IS 

STANT 

IS 

GFTCO 

IS 

OHTCO 

IS 

THE  FUNCTION  TO  BE  MINIMIZED 

ABSOLUTE  ROUGHNESS  OF  TUBING  IN  INCHES 

TUBING  DIAMETER  IN  INCHES 

MONITORED  INTEGRATION  STEP  SIZE  IN  FEET 

GEOTHERMAL  GRADIENT  ,  DEG.  F/  FOOT 

PRESSURE  IN  PSIA 

TEMPERATURE  IN  DEG.  RANKIN 

DENSITY  OF  NAT.  GAS  IN  GM/C.C. 

REYNOLDS  NUMBER 
PRANDTL  NUMBER 
STANTON  NUMBER 

GAS-FILM  HEAT-TRANSFER  COEFFICIENT 
THE  OVERALL  HEAT-TRANSFER  COEFFICIENT 


IMPLICIT  REAL  *8  (A-H,0-Z) 

DIMENSION  PRY  (5)  ,  TR  Y  (5)  ,  PRK  (5)  ,  TRK  (5) 

COMMON  /BL3/AVEMW/BL1 3/DIAM/BL1 4/SG 

COMMON  /BL1 2/PTEMC , ZTEMC , PPREC , Z  PREC 

COMMON  /BL1 5/Q,TSD,PSD, TDD, DEPTH , TSOD, DL 

COMMON  /BL6/CONST/BL 1  8/M 

COMMON  /BL30/VMASS 

COMMON  /BL5C/MMM 

D1=0. 1D-0 1 

D2=0. 1D-C 1 

DL=C .ID  02 

E=C. 6D-03 

RSTAR=X 

RR=E/DIAM 

TGECG=  (TDD-TSOD) /DEPTH 


THE  RUNGE-KUTTA-MERSON  INTEGRATION  ALGORITHM 


TEM=TSD+0 .45967D  03 

PRE=PSD 

T  S  =  T  SOD 

TR=TEM/PTEMC 

PR=FRE/PP  REC 

TRN=TR 

PRN=PR 

TL=C. CD  CO 

WRITE (6,211)  TL,TS,TSD,PRE,ESTAR 

1  DC  14  J 1  =  1 , 5 

IF  ( J 1 - 1 )  2,2,3 

2  PRY  (1)  =PRN 
TRY  (1)  =TRN 
TS=TSOD+TGEOG*TL 

3  IF  (J1-2)  11,4,5 


' 


.  I 


1 
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4  PEY  (2)  =  PR  N  +  PRK  (1) 

TRY  (2)  =TRN+TRK  (1) 

T S=T SOD+TGEOG *  (TL+  DL/0. 3D  01) 

5  IF  (Jl-3)  1%6,7 

6  PRY  (3)  =PEN+C.  5D  0C *PRK  ( 1 ) +0 . 5D  0C*PRK(2) 

TRY (3) =TRN+Q. 5D  0G*TRK  (1 ) +0 • 5D  C0*TRK(2) 
TS=TSOD*TGEOG* (TL+DL/G. 3D  01) 

7  IF  ( J 1 -4)  11,8,9 

8  PRY  (4)  =PRN+0. 375D  00*PRK  (1 )  +0 . 1 1  25D  01*PRK(3) 
TRY  (4) =TRN+0. 375D  00*TRK ( 1 ) +0 . 1 1 25D  01*TRK(3) 

T S=T SOD+TGEOG*  (TL+DL/0. 2D  01) 

9  IF  («J1  -5)  11,10,14 

10  ABC1=PRN+0. 15D  01*PRK(1) 

ABC2=-0.45D  0 1*PRK  (3) +0 . 6D  01*PRK(4) 

PEY  (5) =ABC1+ABC2 

ABC1=TRN+0. 15D  01*TRK(1) 

ABC2=-0.45D  0 1*TRK  (3) +0 • 6D  01*TRK(4) 

TRY  (5)  =ABC 1+ABC2 
TS=TSOD+TGEOG*  (TL+DL) 

PEJ AM  =  PEY  (5) 

TRJAM=TRY  (5) 

11  P  R=PEY ( J 1 ) 

TE=TRY (J1 ) 

TRS=  (TS+0 . 45967D  03) /PTEMC 

PEE=PR*PPREC 

TEM=TR*PTEMC 

C  CALCULATING  IDEAL  HEAT  CAPACITY  OF  NAT,  GAS 

CALL  KASEM  (T R , TEM , SG , M , CPOM) 

TRC=TEM/ZTEMC 

PEC=PRE/ZPREC 

C  CALCULATING  COMPRESSIBILITY  FACTOR, ITS  PARTIAL 

C  DERIVATIVES  W.R.T.  PR  &  TR  , REDUCED 

C  DENSITY,  AND  HEAT  CAPACITY  DEPARTURE 

C  OF  NAT.  GAS 

CALL  BWR  (TEC , PRC , ITER , Z , DR , CPDM , DZDTR , DZDPR) 
IF  (ITER)  22,22,12 

12  DENS=0 • 1 4 9256 D-02*A VEMW*PRE/ (Z*TEM) 
CPM=CPDM+CPOM 

C  CALCULATING  VISCOSITY  OF  NAT.  GAS  IN  C.P. 

CALL  LEE  (TEM, AVEMW, DENS, VISM) 

REN  =  0. 200  1485 5D  02*Q*SG/  ( VISM*DI AM) 

C  CALCULATING  THE  HYDRAULIC  FRICTION  FACTOR 


CALL  MOODY  (E R , REN , FF , I) 


■ 

■ 

' 


on  on 
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C  CALCULATING  THERMAL  COND.  OF  NAT.  GAS 

CALL  MISIC  (TEM,DR, CONST, THKM) 

C  CALCULATING  THE  PRANDTL  NUMBER 

PRAND=CPM*  (VISM/0. ID  03) /  (AVEMW*THKM) 

CALCULATING  THE  STANTON  NUMBER  AND  THE 
GAS-FILM  HEAT-TRANSFER  COEFFICIENT 

B1=FF/0.8D  01 
B2=DSQRT (El) 

ABC3=DSQRT  (PRAND) 

ABC4  = (B1*REN*REN) **G. 27D  CO 
PSI=ABC3*  (0. 1 55D  00*ABC4+C. ID  01/B2) 
STANT=B2/ (PSI-C.45D  01) 
GFTCO=VMASS*CPM*STANT 
U 1 =0 . ID  0 1/GFTCO+RSTAR 
OHTCO=0.1D  01/U1 

SOLVING  FOR  THE  RATE  OF  CHANGE  OF  REDUCED 
PRESSURE  AND  TEMPERATURE  WITH  LENGTH 

CALL  QUON  (TEM,PRE,TR,PR,TRS, Z , DZ DPR , DZDTR , 
1FF,CPM, OHTCO , VMASS, DPRDL , DTRDL, NASER) 

I F  ( NASER- 1 )  22,13,13 

13  PRK  («J1)  =DL*DPRDL/0.  3D  01 
TRK (J1) =DL*DTRDL/G. 3D  01 

14  CONTINUE 

C  MONITORING  THE  INTEGRATION  STEP  SIZE 

AEC5=PRK (1) -0.45D  Q1*PRK(3) 

ABC6=0.4D  01*PRK  (4) -0.5D  CC*PRK(5) 
EPR=ABC5+ABC6 

ABC5=TRK (1) -0.45D  01*TRK(3) 

ABC6=0.4D  C1*TRK  (4) -0.5D  COMTEK  (5) 

ETR  =  ABC  5+  ABC6 
ETEM  =  DABS  (ETR*PTEMC) 

EPRE=DABS  (EPR*PPREC) 

I F  (EPRE-D  1 )  15,15,16 

15  I F  (ETEM-D2)  17,17,16 

16  DL=DL/0.2D  01 
GO  TO  1 

17  PR=PRJAM-EPR 
TR=TRJAM- ETR 
TL=TL+DL 

IF  (EPRE-D  1/0- 1 6D  02) 

18  IF  (ETEM-D2/0.  16D  02) 

IS  DL  =  0. 2D  0  1  *DL 
2C  PRN=PR 


18,20,20 

19,20,20 


TRN=TR 

PRE=PRN*PPREC 

TEM=TRN*PTEMC 

TED=TEM-C . 45967D  03 

WRITE  (6,211)  TL,TS,TED,PRE, RSTAR 

DEL=DEPTH-TL 

IF  (DABS (DEL) -0. ID-04)  23,23,25 
25  IF  (DEL-DL)  21,1,1 

21  DI=DEL 
GO  TO  1 

22  WRI TE  (6 , 2  27) 

23  TED=TEM-0. 45967D  03 
TEMTL=TED 
F=TED-TEMTL 

211  FORMAT (10X,5F10. 2) 

227  FORMAT  (//, 1QX , '***  THE  COURSE  OF  ITERATIONS  IS’ 
*,'  INTERRUPTED  *** * ,//, 7X , * ***  THE  PROGRAM  IS1 
*,'  GOING  TO  ASSUME  ANOTHER  "RSTAR'’***') 

RETURN 

END 


' 

, 

■ 

. 

. 
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SUBROUTINE  AMEIN 

C****  THIS  SUBROUTINE  CAN  BE  USED  WITH 
C****  MAIN  MODEL  AA  ONLY 


C  THIS  SUBROUTINE  CALCULATES  THE  MOLAL  AVERAGE 

C  PROPERTIES, THE  WICHERT  S  AZIZ  CORRECTED 

C  PSEUDOCRITICAL  PROPERTIES  OF  NAT.  GAS  , 

C  AND  OTHER  PROPERTIES 


C  THE  FOLLOWING  TERMS  BELONG  TO  PURE  GASES 


C 

C 

C 

C 

C 


S  IS  THE  SUTHERLAND  CONSTANT  IN  DEG.  K 

VC  '  IS  CRITICAL  MOLAL  VOLUME,  C.C./GM  MOLE 
V CLC  IS  CRITICAL  VOLUME  ,  CU.  FT/  LB. 

ZC  IS  CRITICAL  COMPRESSIBILITY  FACTOR 
GAMA  IS  CONSTANT  DEFINED  BY  EQ.  58  IN  TEXT 


C  THE  FOLLOWING  TERMS  BELONG  TO  NAT.  GAS 


C 

C 

C 

c 

c 

c 

c 

c 

c 


AVEMW  IS 
SG  IS 
PPREC  IS 
PPC  IS 
PTEMC  IS 
PTC  IS 
GAM  IS 
ZTEMC  IS 
ZPREC  IS 


AVERAGE  MOLECULAR  WEIGHT 
GAS  GRAVITY 

PSEUDOCRITICCAL  PRESSURE  IN  PSIA 
PSEUDOCRITICCAL  PRESSURE  IN  ATMOSPHERS 
PSEUDOCRITICAL  TEMPERATURE , DEG .  R 
PSEUDOCRITICAL  TEMPERATURE , DEG .  K 
CONSTANT  DEFINED  BY  EQ.  58  IN  TEXT 
CORRECTED  PSEUDOCRITICAL  TEMPERATURE 
CORRECTED  PSEUDOCRITICAL  PRESSURE 


IMPLICIT  EEAL*8  (A-H,0-Z) 

DIMENSION  PC  (15) 

COMMON  /BL3/A VEM W/BL4/SEGMA  (15)  ,EOK(15) 

COMMON  /BL5/YMW  (15)  ,TC  (1  5)  ,  ZC  (1  5)  ,G  AM  A  (15) 

CCMMON  /BL9/YMOLF (15) /BLIO/TNBOI (15)  /BL20/S (15) 

COMMON  /BL11/TEMC  (15)  ,PREC  (15)  , VOLC (15)  /BL6/CON ST 

CCMMON  /BL12/PTEMC, ZTEMC, PPREC, ZPREC 

CCMMON  /BL14/SG/BL18/M 

CCMMON  /BL50/MMM 

COST 1 =0 • 1 5D  01/0. 18D  01 

CCST2=0 . ID  01/C.6D  01 

COST3=0.2D  01/0. 3D  01 

Y 1 1 = YMOLF  (1) 

Y22=YMOLF  (2) 

Y33  =YMOLF  (3) 

IF(M-I)  3,1,1 
1  AVEMW=0.0D  00 
PTEMC=0.0D  00 
PPREC=0 . 0 D  00 
DO  2  1=1,15 
S  (I)  =COST1*TNBOI  (I) 

PC  (I) =PREC  (I) /C.  146964757D  02 


■ 


. 


' 


■ 


. 
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TC  (I)  =  TEMC  (I)  /O.  18D  01 

VC=0. 6238 268 96D  02*YMW  (I) *VOLC  (I) 

ZC  (I)  =PC  (I)  *  VC/  (0*  8206D  02*TC(I)  ) 

Y  1  =  TC (I) *  *C0ST2 

Y  2  =  PC (I) * *C0ST3 

GAMA  (I) =0 . ID  07*Y1*DSQRT (YMW (I)  )  /Y2 
A  VEMW=A  VEMW+ YMOLF (I) *YMW (I) 

PTEMC=PTE  MC+ YMOLF (I) *TEMC (I) 

2  PPREC=PPR EC+ Y MOLF  (I)  *PREC  (I) 

IF(MMM-I)  10,3,3 

10  SG=A VEM W/C .  2897D  02 
GO  TO  7 

3  AVEMW=SG*0.2897D  02 

C  THE  BROWN  ST  AL  AND  CARR  ET  AL  CORRELATIONS 

C  FOP.  THE  PSEUDOCRITICAL  PROPERTIES  OF 

C  SWEET  AND  SOUR  NAT.  GASES 

PTEMC=0.171137D  03+0.313725D  03*SG 
TCCOR=-0. 8D  0  2*Y1 1 -0 .  25D  03*Y22+0.13D  Q3*Y33 
PTE  MC=PTEMC+TCCOR 
IF (SG-C.85D  00)  4,4,5 

4  PPREC=C • 6951D  03-0. 4D  02*SG 
GO  TO  6 

5  PPREC=C . 7 043 96D  C3-0.  51724D  02*SG 

6  PCCOR=0 . 4  4 D  03*Y 1 1-0. 17D  03*Y22+0.6D  03*Y33 
PPREC=PPREC-*-PCCOR 

C  THE  WICHERT  &  AZIZ  CORRECTED  PSEUDOCRITICALS 

C  DUE  TO  PRESENCE  OF  C02  AND  H2S 

7  S 1  =  YMOLF ( 1 ) ♦ YMOLF  (3) 

S2  =  YMOLF  (3) 

E3 1 =0. 1  2D  03* (SI **0 . 9D  00-S1**0.  16D  01) 
E32=0.15D  02*  (DSQRT  ( S 2) -S2**4) 

E3=E31+E3  2 

S4=PTEMC+S2* (0. ID  01-S2)*E3 

ZTEMC=PTEMC-E3 

ZPREC=PPREC*ZTEMC/S4 

PTC=PTEMC/0. 1 8D  01 

PPC=PPREC/0. 14696476D  02 

G AM=PTC**COST2*DSQRT (AVEMW) /PPC**COST3 

CCNST=GAM*0. 14348907D  06 

RETURN 

END 


■ 


SUBROUTINE  HISIC  (TEM  , DR , CONST, THKM) 

C****  THIS  SUBROUTINE  CAN  BE  USED  WITH 
C****  MAIN  MODEL  AA  ONLY 


C  THIS  SUBROUTINE  CALCULATES  THE  THERMAL 

C  CONDUCTIVITY  OF  GAS  MIXTURES 

C**  THE  HISIC  AND  THODOS  METHODS  ARE  USED  TO 
C  CALCULATE  THKO,THKMO, THKMD, AND  THKM 

C  TEM  IS  TEMP.  IN  DEGREE  RANKIN 

C  DR  IS  REDUCED  DENSITY  OF  GAS  MIXTURE 

C  THKMO  IS  THERMAL  COND.  OF  GAS  MIX.  AT  C  PRESS. 

C  THKMD  IS  THERMAL  COND.  DEPARTURE 

C  THKM  IS  THERMAL  CONDUCTIVITY  OF  GAS  MIXTURE 

C  IN  CAL/CM  SEC  K 

C  CONST  IS  1C**8*0. 27**5*GAM 


C 

C 

C 

C 

C 

C 

C 

C 


THE  FOLLOWING  TERMS  BELONG  TO  PURE  GASES 
Y MW  IS  MOLECULAR  WEIGHT 

TC  IS  CRITICAL  TEMPERATURE  IN  DEG.  KELVIN 
ZC  IS  CRITICAL  COMPRESSIBILITY  FACTOR 
S  IS  SUTHERLAND’S  CONSTANT  IN  DEG.  KELVIN 

CPO  IS  HEAT  CAPACITY  AT  ZERO  PRESSURE 
IN  CAL/GM  MOLE  K 

THKO  IS  THERMAL  COND.  AT  0  PRESSURE 
YMCLF  IS  MOLE  FRACTION 


IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  A  (15,15)  ,SS  (15, 15)  ,VISO  (15)  , 

*  THKC(15)  ,F(2)  ,TRED(15) 

COMMON  /BL4/SEGMA  (15)  ,EOK  (15) 

COMMON  /BL5/YMW  (1  5)  ,  TC  (1  5)  ,  ZC  (1  5)  , GAMA  (15) 

COMMON  /BL7/CPO (15) /BL9/YMOLF (1 5) /BL20/S ( 1 5) 

UNITY=0.1D  01 

FRAC=2.0D  00/3. 0D  00 

TK=TEM/0. 18D  01 

TKK=TK 

DO  8  1=1 ,  1  5 

TRED  (I)  =T K/TC  (I) 

DO  5  K  =  1 , 2 
IFLAG=K 
T=TK/EOK (I) 

IF  (T-0 • 2D  01)  1,2,2 

1  F  (K) =0. 634D  CO*  (UNITY+0.475D  C0*DLOG(T)) 

GO  TO  3 

2  F (K) =0 •  69 7D  00*  (UNITY+C.  323D  CO*DLOG(T)) 

3  IF  (TRED  (I) -0. 2D  01)  4,6,6 

4  TK=C . 2D  0  1*TC  (I) 

5  CONTINUE 

Y  =  DSQRT  (YMW (I) *TK) 


■ 


A  1 


6  VISO  (I) =0. 26693D-02*Y*F (IFLAG) /SEGM A (I)  **2 
IF  (IFLAG- 1 )  8,8,7 

7  TK=TKK 

VISO  (I)  =VISO  (I)  *0. 5D  00*TRED  (I)  *F  (1 )  /F  (2) 

8  CONTINUE 

DO  11  1=1,15 
IF  (1*3)  9,9,10 

9  A  1=0 . 1 9 5D  G3*ZC  (I) -C. 3194D  02 
A 2=0 ,16  83d  02-0. 825D  02*ZC (I) 

A3=0. 1  524  D  01-0.  28D  C1*ZC(I) 

A 4=  (A1*TRED (I) +  A2) **A3 

THKO  (I) =CPO  (I) **0. 75D  CQ*A4/GAMA  (I) 

GO  TO  11 

10  B  1=0. 1 452D  02*TRED  (I) -0. 514D  01 
THKO  (I)  =CFO  (I)  *B1**FRAC/GAMA  (I) 

11  CONTINUE 

DO  12  1=7,15 
DC  12  J=1 ,15 
SS (I,  J)  =D SQRT (S (I) *S  (J) ) 

C 1  =  (UNITY  +  SS  (I,  J)/TK)  /  (UNITY* S  (I)  /TK) 

C2=  (VISO  (I) /VISO  (J)  )  *  (YMW  (J)/YMW  (I)  )  **C.75D  00 
C 3=  (UNITY  +  S  (I) /TK)  /(UNITY+S  (J)/TK) 

A  (I,J) =0. 25D  G0*C1*  (U  NITY  +  DSQRT ( C2*C3 ) ) **2 

12  CONTINUE 
THKMO=G.OD  00 
DO  15  1=1 ,15 

IF  (YMOLF (I) -0. ID-06)  15,15,25 

25  DUM Y=G • ID  01 
DC  14  j=1 ,15 

IF  (YMOLF  (J) -0. 1D-C6)  14,14,26 

26  IF(I-J)  13,14,13 

13  DUM Y=DUMY+ A  (I, J) * YMOLF (J) /YMOLF (I) 

14  CONTINUE 

THKMO=THK  MO+THKO (I)  /DUMY 

15  CONTINUE 

IF  (DR-0.5D  00)  16,16,17 

16  X1=DEXP (0.535D  00*DR) 

T HKMD=0 . 1 4D  02*  (X1-UNITY) /CONST 
GO  TO  20 

17  IF  (DR-C .2D  01)  18,18,19 

18  X  1  =  DEXP  (0 • 67D  CQ*DR) 

THKMD=C . 1 3  ID  02*  (X 1 -0 . 1 069D  01)/CONST 
GO  TO  20 

19  X1=DEXP  (0. 1155D  01*DR) 

THKMD=C. 2976D  0 1  *  (X 1  +  0. 20 1 6D  01) /CONST 
2C  THKM=THKMO+TH KMD 
RETURN 
END 


■ 

■ 

■ 

■ 

■ 

■  • 

' 

■ 

. 

. 
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S UEROUTINE  SEARCH  (F ,EPS, X, YMIN, YMAX , RSTOUT) 

C****  THIS  SUBPOOTINE  CAN  BE  USED  WITH 
C****  EITHER  MAIN  MODEL  AA  OR  BB 

C  THIS  SUBROUTINE  SEARCHES  FOR  THE  VALUE  OF 

C  B*  OR  R**  WHICH  MINIMIZES  THE  FUNCTION  'F' 

IMPLICIT  REAL*8  (A-H,0-Z) 

EXTERNAL  F 

DIMENSION  RST  (3)  , TEST  (2)  ,DELTA  (3) 

RST  (1)  =  YMIN 
RST  (3)  =  Y M A X 
ML=0 

30  IF  (ML- 1 7)  1,1,21 

21  W RITE  (6,2  25) 

GO  TO  13 

1  RST  (2)  =  (RST  (1) +RST  (3)  )/0.  2D  01 

IF (RST (2) -RSTOUT)  13,13,2 

2  DO  6  K3=1 ,3 
IF  (ML)  4,4,3 

3  IF  (K3-2)  6,4,6 

4  RSTAR=RST ( K 3 ) 

X=RSTAR 
DIFF=F (X) 

IF (DABS (D IFF) -EPS)  13,13,5 

5  DELTA  (K3)  =DIFF 

6  CONTINUE 

TEST  (1)  =DELT A  (1)  *DELTA  (2) 

TEST  (2)  =DELTA  (3)  *DELTA  (2) 

DESCN=DELTA (1) *DELTA  (3) 

IF  (ML)  23,23,24 

23  IE(DESCN)  24,  24,25 
25  WRITE (6,226) 

GO  TO  13 

24  ML=ML+ 1 

DO  8  K 5  =  1 ,2 
IF  (TEST  (K 5 ) )  7,13,8 

7  ITEST=K5 
GC  TO  9 

8  CONTINUE 

9  I F (ITEST-  1 )  10,10,11 

10  RST  (3)  =RST  (2) 

DELTA (3 ) =  DELTA  (2) 

GO  TO  30 

11  IF  (TEST  (2))  12,13,12 

12  RST  (1)  =RST  (2) 

DELTA  (1)  =DELTA  (2) 

GO  TO  30 

225  FORMAT  (// ,20X,****  NO  SOLUTION  IS  ', 

*  *  REACHED  WITH  20  ITERATIONS  ***') 

226  FORMAT  (//,20X,'***  NO  SOLUTION  IS  ', 


4  • 

' 


*  *  PROVIDED  TO  THIS  PROBLEM  ***') 
13  RETURN 


END 

I 


fN  no  cr 
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SUBROUTINE  KASEM  (TB , TEM , SG , M , CPOM) 

C****  THIS  SUBROUTINE  CAN  BE  USED  WITH 
C****  EITHER  MAIN  MODEL  AA  OR  BB 


c 

TR 

IS 

REDUCED  TEMP. 

c 

TEM 

IS 

TEMPERATURE  IN  DEGREE  RANKIN 

c 

SG 

IS 

GAS  GRAVITY 

c 

M 

IS 

A  FLAG  INTEGER 

c 

M=0 

:  SUPPLY  SG,H2S,C02,AND  N2 

c 

M  =  1 

:  SUPPLY  GAS  ANALYSIS 

c 

CPOM 

IS 

HEAT  CAPACITY  AT  ZERO  PRESS. 

OF 

NAT. 

GAS 

c 

CPC 

IS 

HEAT  CAPACITY  AT  ZERO  PRESS. 

OF 

PURE 

GAS 

c 

YMOLF 

IS 

MOLE  FRACTION 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  A  (1 4) 

COMMON  /BL7/CPO  (15) /BL8/CPOC (15,4) /BL9/YMOLF (15) 
DATA  A/0. 5596695D  0 1 , -0 . 22334 80 D  01, 

1  0. 807265D  €0,-0. 10C390D  01,0.31416CD  01, 

0. 5758700D  01,0.425860D  00 ,0 . 1 24323D-02 , 
-0.4C560D-01 ,0. 100889D-02,0. 3623622D  00, 
-C.4661581D  00  ,0. 97  57D-0 3 ,0. 2708  21 67D-Q2/ 
IF(M-I)  4,1,1 

1  CPOM=0. CD  00 
DO  3  1  =  1 ,  5 
CPO  (I) =0. CD  00 
DO  2  J=1, 4 

2  CPO  (I)  =CPO  (I)  +CPOC  (I,  J)  *TEM**  (J-1) 

3  CRCM=CPOM+CPO  (I) *YMOLF (I) 

GC  TO  5 

4  B=A  (1)  +  A  (2)*SG  +  A  (3)  *SG**2 
C  =  A  (4)  +  A  (5)  *SG+A  (6)  *SG**2 
CPCM=B+C*TR 

FCO2=0 .  ID  Q1+YMOLF  (1)  *  (A  (7)  +A  (8)  *TEM) 

FH2 S=0 • ID  OI  +  YMOLF  (3)  *(A  (9) +A  (10) *TEM) 

Y  1  =A  (11)  * YMCLF  (2)  +  A  (12)  * YMOLF  (2)  **2 
Y 2  =  A  (13)  * YMOLF  (2)  +A  (14)  *YMOLF  (2)  **2 
FN2  =  0. ID  0 1  + Y 1  +Y2*TEM 
CPOM=CPOM/  (FC02*FN2*FH2S) 

5  RETURN 
END 


' 


SUBROUTINE  BWR  (TR , PR , ITER , Z , DR ,C PDM , DZDTR , DZDPE ) 


C****  THIS  SUBROUTINE  CAN  BE  USED  WITH 
C****  EITHER  MAIN  MODEL  AA  OR  BB 


W.R.T.  TR 

W.R.T.  PR 


THE  PSEUDO-REDUCED  TEMPERATURE 
THE  PSEUDO-REDUCED  PRESSURE 
THE  GAS  COMPRESSIBILITY  FACTOR 
THE  REDUCED  DENSITY  OF  GAS 
THE  PARTIAL  DERIVATIVE  OF  Z 
THE  PARTIAL  DERIVATIVE  OF  Z 
THE  HEAT  CAPACITY  DEPARTURE 
AN  ITERATION  FLAG 
EITHER  OR  BOTH  TR  AND  PR  ARE 
OUTSIDE  RANGE  OF  CORRELATION 
NUMBER  OF  ITERATIONS 


REAL*8  (A- H, O-Z) 

DIMENSION  A  (8) 

DATA  A/C.  31506237D  CO  ,  1 . 0467G99D  00, 

A  0. 57832729D  00 , 0 . 5 3530771 D  00 ,0.  6 1 232032D  00, 
B  0. 10488813D  00 , 0 . 6 81 5700 1 D  00 , 0. 68446 549D  00/ 
ITEB=0 
DR=C‘ •  ID  01 


c 

TR 

IS 

c 

PR 

IS 

c 

Z 

IS 

c 

DR 

IS 

c 

DZDTR 

IS 

c 

DZDPR 

IS 

c 

CPDM 

IS 

c 

ITER 

IS 

c 

ITER 

=0 

c 

c 

ITER>0 

IMPLICIT 

IF  (TR-0. 1C5D  01)  10, 1 ,  1 

1  IF (TR-0.3D  01) 2,2,10 

2  IF  (PB-0 . 30D  02)3,3,10 

3  DO  8  ITER= 1 , 30 
DR2=DR*DR 

T1=  (A  (1) *TR- A  (2) -A (3) /TR**2) *DR 
T2=  (A  (4)  *TR-A  (5)  )  *DR2 
T3  =  A  (5)  *A  (6)  *DR**5 
T4  =  A  (7) *DR2/TR**2 
T 5  =  A  (8)  *DR2 
T6  =  DEXP  (-T5) 

P=  (TR+T1+T2+T3) *DR+T4*DR*  (0. ID  01+T5) *T6 
DPI =TR+G. 2D  0 1 *T1 +0*  3D  G1*T2+0.6D  01*T3 
DP2=T4*T6*  (0. 3D  01+0. 3D  C1*T5-0.2D  01*T5*T5) 
DP  =  EP1  +  DP  2 

DR1=DR- (P-0.270D  00*PE)/DP 
IF  (DR 1 ) 4,4,5 

4  DR1 =0 . 5  D  C  0 * D R 

5  IF (DR1-0.22D  01) 7,7,6 

6  DR1 =DR+0. 9D  00*(Q.22D  01-DR) 

7  IE (DABS  (DR-DR1) -0. ID-04) 9,8, 8 

8  DR=DR1 

9  Z  =  0 • 270D  C 0*PR/  (DR  1 *TR) 

T7=  (A  (1) *TR+0. 2D  G1*A  (3) /TR**2) *DR 
T 8  =  A  (4) *TR*DR**2 
T9  =  A  (3)  *DP/TR**2 
T10=A  (7)  /A  (8)  /TR**2 
R 1 =DP/TR 


■ 


R2=  (TR+T7+T8-0.2D  £}  1  *T4*T6*  (0 • 1 D  01+T5))/TR 

F2=B2*DR 

P3  =  0. 270D  00/DP 

Y 1 =0 . 3D  Q1*T6*(0.2D  01*T10+T4) 

R3=  (-TR+9.6D  01*  (T9-T10) +Y1) /TR 
P4=P2/DP 

DZDTR=Z* (P4/DR-0. ID  01/TR) 

DZDPF.  =  Z*  (0.  ID  01/PR-P3/DR) 

CPDM=0.  19 8654284D  01*  (R3  +  R2**2/R 1 ) 

10  RETURN 
END 
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SUBROUTINE  LEE  (TEM , A VEM W , DE NS , VI SM) 

C****  THIS  SUBROUTINE  CAN  BE  USED  WITH 
C****  EITHER  MAIN  MODEL  AA  OR  BB 

C  THIS  SUBROUTINE  IS  THE  LEE, ET  AL.  VISCOCITY  EQ. 

C  TEM  IS  TEMPERATURE  IN  DEGREE  RANKIN 

C  A VEMW  IS  AVERAGE  MOLECULAR  WEIGHT 

C  DENS  IS  THE  GAS  DENSITY  IN  GRAM/C. C. 

C  VISM  IS  THE  GAS  VISCOCITY  IN  C. P. 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  V(11) 

DATA  V/0.94D  01 ,0. 2D-01 ,0. 1 5D  01,0.209D  03, 

1  Q.19D  02,0. 35D  01,0.986D  03,0. ID-01, 

2  0.24D  01,0. 2D  00,0. ID  05/ 

VKl  =  (V  (1)  +V  (2)  *AVEMW)  *TEM**V  (3) 

VK2  =  V  (4)  +  V  (5)  *AVEM W+TEM 

VK= VK1/VK2 

VX  =  V  (6)  +V  (7)  /TEM  +  V  (8)  *AVEMW 
V  Y= V (9) -V (10) *VX 
VY1=VX*DENS**VY 
VISM=VK*DEXP  (VY1) /V  (1 1) 

RETURN 

END 
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SUBROUTINE  MOODY  (RR , REN, FF, I) 

C****  THIS  SUBROUTINE  CAN  BE  USED  WITH 
C****  EITHER  MAIN  MODEL  AA  OR  BB 

C  RR  IS  THE  RELATIVE  ROUGHNESS 

C  REN  IS  REYNOLDS  NUMBER 

C  RENCR  IS  THE  BOUNDARY  BETWEEN  THE 

C  TRANSITION  AND  ROUGH  PIPE  REGIONS 

C  FF  IS  THE  HYDRAULIC  FRICTION  FACTOR 

C  I  IS  NO.  OF  ITERATIONS  REQUIRED 

C  1=0  ROUGH  PIPE  REGION 

C  I>0  TRANSITION  REGION  . 

IMPLICIT  PEAL *8  (A-H,0-Z) 

DIMENSION  G  (1 5) 

1  =  0 

C  =  DLOG  (0.  ID  02) 

IF  (RR)  2,2,1 

1  RENCR=G . 3  2D  04/RR 
IF  (REN-RENCR)  2,6,6 

C  THE  COLEBROOK  TRANSITION  LAW 

2  X  =G . 8D-C  1 
DO  4  1=1, 15 

X 1  =  RR/0 . 3  7D  01+0. 251D  01/(REN*X) 

X2  =  DLOG 1C  (XI) 

FX=0.1D  01+0. 2D  C1*X*X2 

FXD  =  0. 2D  0 1 *X  2-0 . 50  2D  0 1/ (REN*X*X1*C) 

G  (I) =X-FX/FXD 
X=G  (I) 

IF(I-I)  4,4,3 

3  IF  (DABS  (G  (I) -G  (1-1)  ) -0.  ID-05)  5,5,4 

4  CONTINUE 

5  FF=X*X 
GO  TO  7 

C  VON-KARMAN'S  ROUGH  PIPE  FORMULA 

6  FF  =  C • 1 D  01/(0. 2D  01 *DLOG1C  (0. 37D  01/RR) ) **2 

7  CONTINUE 
RETURN 
END 


. 


non 
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SUBROUTINE  QUON (TEM , P EE , TR , PE ,TR S, Z , DZD PE , 
1DZDTR,FF, CPM,OHTCO,VM ASS, DPRDL, DTRDL, NASER) 

c****  this  SUBFOOTINE  CAN  BE  USED  WITH 
C****  EITHER  MAIN  MODEL  AA  OR  BB 

THIS  SUBROUTINE  SOLVES  FOR  THE  RATE  OF  CHANGE 
OF  REDUCED  PRESSURE  &  TEMP.  WITH  LENGTH 
(SEE  EQUATIONS  6  &  7  IN  TEXT  ) 

IMPLICIT  REAL*8  (A-H,0-Z) 

COMMON  /BL1/GASUC ,GC , AJ, S LOPE /BL 2/1 TYPE 
COMMON  /BL3/A VEM W/BL1 3/DIAM 
N  ASER=1 
0=SL0PE 
BETA=0.1D  01 
DIAN=DIAM/0. 12D  02 
A1=EETA*VMASS**2*TEM*GASUC 
PEF=PRE*Q • 14 4 D  03 
A2=GC*AJ*PRF**2 
A  3=A 1/A2 
A<i  =  A3*TEM*GASUC 
A5=Z/PR-DZDPR 
A6=Z/TR+DZDTR 

ASTAR=-A1 *A5*PRF+GC*PRF**3/  (PR*AVEMW) 
BSTAR=A1*A6*PRF 

CSTA1=-A2*PRF*0/  (Z*TEM*AJ*GASUC) 

CSTA2=-A1 *Z*PRF*FF/  (0.2D  0 1 *DIAN*BETA) 

C  ST AR=C  STA1+CSTA2 

DSTAR»A4*A5*Z+TEM*TR*GASUC*DZDTR/(AJ*PR*AVEMW) 
ESTAR=-A4*A6*Z-TEM*CPM/  (TR*AVEMW) 
CHRIS=VMASS*DIAN 
FSTA1=CHRIS*0 

F  STA2  =  0 . 4 D  0 1 *OHTCO*TEM*A J*  (TRS-TR) / (TR* A VEM W) 
FST AR=  (FSTA1 -FSTA2) /  (CHRIS* A J) 
A7=ASTAR*ESTAR-BSTAR*DSTAR 
IF  ( A7)  1,3,1 

1  DPEDL=  (CSTAR*  ESTAR-BSTAR*FSTAR) /A 7 
DTRDL=  (ASTAR*FSTAR-CSTAR*DSTAR) /A7 
IF(ITYPE-I)  4,2,2 

2  DPRDL=-DP RDL 
DTRDL=-DTEDL 
GC  TO  4 

3  N ASER=0 

4  RETURN 
END 


' 

. 


C*******************  MAIN  MODEL  BB  **************** 
C*************************************************** 
C**  THIS  MAIN  MODEL, 'BB',  REQUIRES  THE  FOLLOWING 
C**  FUNCTION  AND  SUBROUTINE  PROGRAMS  : 

C**  F  &  AMEIN  (LISTED  IN  PAGES  TO  FOLLOW)  , 

C**  SEARCH, KASIM, BWR, LEE, MOODY  &  QUON  (LISTED 

C**  IN  PREVIOUS  PAGES)  . 

C**  HOW  TO  USE  THIS  COMPUTER  PROGRAM 


C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 


ITYPS  IS  A  FLAG  INTEGER 

IT YPE=G  :  INJECTION  PROBLEM 
ITY PE= 1  :  PRODUCTION  PROBLEM 
M  IS  A  FLAG  INTEGER 

M=0  :  SUFPLY  SG  ,H2S ,C02 ,AND  N2 
M=1  :  SUPPLY  GAS  ANALYSIS 

MM  IS  A  FLAG  INTEGER 

MM=C  :  DO  NOT  SUPPLY  PS EUDOCRITICALS 
MM=  1  :  SUPPLY  PSEUDOCRITICALS 
MMM  IS  A  FLAG  INTEGER 

MMM=0  :  ONE-PHASE  FLOW 
MMM=1  :  TWO-PHASE  FLOW 


C  COMBINATIONS  OF  THE  FLAG  INTEGERS  ARE  AS  FOLLOWS: 


C 

ITYPE 

M 

MM 

MMM 

DESCRIPTION 

C 

1 

1 

1 

1 

TWO-PHASE  FLOW  PRODUCTION 

c 

PROBLEM.  SUPPLY  GAS  COMPOSITION 

c 

AND  SG  OF  FLOWING  STREAM 

c 

1 

1 

1 

0 

ONE-PHASE  FLOW  PRODUCTION 

c 

PROBLEM.  SUPPLY  GAS  COMPOSITION 

r* 

u 

1 

0 

0 

1 

TWO-PHASE  FLOW  PRODUCTION 

c 

PROBLEM.  SUPPLY  SG  OF  FLOWING 

c 

STREAM, C02,N2  &  H2S 

c 

1 

0 

A 

0 

sJ 

ONE-PHASE  FLOW  PRODUCTION 

c 

PROBLEM.  SUPPLY  SG,C02,N2  &  H2S 

c 

1 

0 

1 

0 

ONE-PHASE  FLOW  PRODUCTION 

c 

PROBLEM.  SUPPLY  SG , C02 , N2 , H2  S 

c 

AND  THE  REPORTED  PTEMC  &  PPREC 

c 

GASUC 

IS 

GAS 

UNIVERSAL  CONSTANT 

c 

GC 

IS 

CONVERSION  FACTOR 

c 

A  J 

IS 

CONVERSION  FACTOR 

c 

SLOPE 

IS 

SLOPE  OF 

PIPE  SEGMENTS 

C**  THE  FOLLOWING  TERMS  BELONG  TO  PURE  GASES 
C  YMOLF  IS  MOLE  FRACTION 

C  Y MW  IS  MOLECULAR  WEIGHT 

C  YMW1  IS  MOLECULAR  WEIGHT 

C  T EMC  IS  CRITICAL  TEMPERATURE  IN  DEG.  RANKIN 


' 

■ 


. 


Cr  UJ  NJ  —k  C  OJ  M  -*  -CrOJNJ 
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c 

TEMC1 

IS 

c 

PREC 

IS 

c 

PREC1 

IS 

c 

CPOC 

IS 

c 

c 

c 

CPOC1 

IS 

c 

SG 

IS 

c 

PTEMC 

IS 

c 

P  PREC 

IS 

c 

DIAM 

IS 

c 

DEPTH 

IS 

c 

TSOD 

IS 

c 

TSD 

IS 

c 

TDD 

IS 

c 

PSD 

IS 

c 

Q 

IS 

IMPLICIT 

EXTERNAL 

CRITICAL 
CRITICAL 
CRITICAL 
THE  FOUR 


TEMPERATURE  IN  DEG.  RANKIN 
PRESSURE  IN  PSIA 
PRESSURE  IN  PSIA 
COEFFICIENTS  OF  THE  IDEAL 


HEAT  CAPACITY  EQ. 
(KOBE  ET  AL.  ) 

THE  SAME  AS  CPOC 


OF  PURE  GASES 


1 


GAS  GRAVITY 

PSEUDOCRITICAL  TEMP.  IN  DEG.  RANKIN 
FSEUDOCRITIC AL  PRESS.  IN  PSIA 
TUBING  DIAMETER  IN  INCHES 
LENGTH  OF  TUBING  IN  FEET 
GROUND  SURFACE  TEMP.  IN  DEG.  F 
TUBING  HEAD  TEMP.  IN  DEG.  F 
BOTTOM-HOLE  TEMP.  IN  DEG.  F 
TUBING  HEAD  PRESSURE  IN  PSIA 
TOTAL  FLOW  RATE  IN  MSCF/DA Y 

REAL*8  (A-H,0-Z) 

F 

DIMENSION  YMW1  (15)  ,TEMC1  (15)  ,PREC1  (15) 

DIMENSION  CPOC1  (15,4) 

COMMON  /BL1/GASUC,GC, AJ, S LOPE /BL 2/1 TYPE 
COMMON  /BL8/CPOC (1 5, 4)/BL5/YMW (1 5) 

CCMMGN  /BL11/TEMC (15)  ,PREC  (15) 

COMMON  /BL12/PTEMC, ZT EMC, P PREC, ZPREC/BL18/M 
COMMON  /BL 1 3/DIAM/BL1 4/SG/BL9/YMOLF (15) 

COMMON  /BL15/Q,TSD,PSD,TDD,DEPTH,TS0D,DL 
COMMON  /BL30/VMASS/BL31/MM 
COMMON  /B150/MMM 

DATA  YMW  1/0.  44013D  C2,0.28013D  02,0.340760  02, 

0 2, C . 30C70D  02,0. 44097D  02, 
02,0. 5812 4D  02,0.72151D  02, 
02,0. 86178D  02,0.86178D  02, 
03,0. 1CC21D  03, 0. 1 1 423D  03/ 
DATA  TEMC1/0. 54757D  03,G.22727D  03,0.67237D  03, 

03, 0. 66568D  03, 


YMW 1/0.  440 1  3D 
0.  1 604  3D 
0. 58124D 
0. 72151D 
0. 10021 D 
TEMC1/0. 54757D  03,G.22727D 
C.  34304D  0 3 ,0 . 54976D 


0.73465D  03,0. 76532D  03,0.82877D  03, 
0.84537D  03,0. 89550D  03,0.91337D  03, 
0.95467D  03,0. 97247D  03,0.10239D  04/ 
PR EC  1/0 . 1 0 7 1 D  04,0. 4930D  03,0.  1306D  04, 
0.6678D  03 ,0. 7C78D  03,G.6163D  03, 
03,0. 5507D  0  3,0.4 904 D  03, 
03,0. 4366D  03 , 0. 4369D  03, 
03,0.3968D  03,0.3606D  03/ 

01,  0.6903D  01,  0.7070D  01, 
01,  C.1648D  01,-0. 96 60D  00, 
01,  0.9450D  00,-0. 2273D  01, 

0 1 ,-0 • 908 1 D  00,  0.1657D  01, 
01,  0.31899D  01,  0.35844D  01 


DATA 


0. 529 1 D 
0. 4886D 
0. 3965D 
DATA  CPOC1/0. 5316D 
0.  4750  D 
-0.  1890D 
0.  1618D 
-0. 26C69D 


0. 793611 1D-C 2, -0.20 85000D-03,  0. 1737777D-G2 


/ 

/ 
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F 

G 

H 

I 

J 

K 

L 

M 

0 

P 

Q 

R 

S 


0. 6666666D-02, 
0. 5519999D-01 , 
0. 6027777D-01, 
0. 96 1 105  5D-0 1 , 
-0.2580864D-05, 
0. 9351 851 D-06 , 
-0. 1695S87D-04, 
-0. 1655864D-04, 
-0. 3089506D-04, 
0.3058984D-G9, 
-0. 45C9602D-09, 
0 • 20  4389  5D-08 , 
Q. 1731824D-08, 
0. 3834362D-08, 
GC=0 . 32 1 74D  02 
GASUC=0. 1 545D  04 
AJ  =  C. 7781 6D  03 


C  .  2291 1 1 1D- 01 , 
0 • 4929444D-G  1 , 
0.7843494D-Q1 , 
C • 83 1 5605D-Q 1 , 
0 • 5956790D-C6, 
•0  •  4722222D-G  5,- 
■0. 1351851D-04,- 
■0  •  239  69 1 0D-O  4 ,  ■ 
•0  •  237  31  48D-04  ,  ■ 
•0. 11764  40  D- 0  9,- 
0 .  2983539D-09, 
C.  1433470D-08, 
0  •  2848422D-0  8, 
0.261 4881D-0  8, 


0. 4043888D-0 1 , 
0.  6907777D-C  1 , 
0.73  27777D-0  1  , 
0. 9502444D-0  1 , 
0. 4209876D-0  6 , 
0.11 58950D-C4, 
0. 21 90432D-0  4, 
•0.21 12345D-0  4 , 
■0. 2742592D-C4 , 
0. 1348936D-09, 
0. 129972 5D- 0  8, 
0. 2719478D-08, 
0.23  62  825D-0  8 , 
0. 306001 2D-C  8/ 


SIOPE=0.1D  01 
DO  5  1=1 , 15 
Y MW  (I)  =  YM W 1  (I) 

TEMC  (I)  =T EMC  1  (I) 

PEEC  (I) =PREC1  (I) 

DO  5  J=1 , 4 

5  CPOC  (I, J) =CP0C1  (I, J) 

WRITE  (6,208) 

WRITE  (6,2^2)  (YMW  (I)  ,  TEMC  (I)  ,  PREC  (I)  ,1=1,15) 
WRITE  (6,213) 

WRITE  (6,205)  (  (CPOC  (I, J)  , J  =  1 ,4)  ,1  =  1 , 15) 

WRITE  (6,214) 


C  READ  3  CARDS.  GAS  COMPOSITION. EACH  CARD  HAS 

C  THE  MOLE  FRACTIONS  OF  5  COMPONENTS.  THE 

C  ORDER  OF  COMPONENTS  IS  AS  FOLLOWS  : 

C  C02,N2,H2S,C1,C2,C3,I-C4,N-C4,I-C5,N-C5, 

C  I-C6, N-C6 ,I-C7,N-C7,N-C8  . 

C  IF  GAS  COMPOSITION  IS  NOT  AVAILABLE, 

C  SUPPLY  3  BLANK  CARDS  . 


READ  (5, 20  6)  (YMOLF  (I)  ,1=1 , 1 5) 

WRITE  (6,215) 

WRITE(6,2C7)  (YMOLF  (I)  ,1  =  1 , 1 5) 

C  READ  1  CARD  (  COMBINATION  OF  FLAG  INTEGERS 

C  AND  S G , PTEMC  8  PPREC  ) 

READ  (5,20  9)  I TYPE , M , MM , MM M , SG , PTEMC , PPREC 

C  READ  2  CARDS  (  WELL  PRODUCTION  DATA  ) 

READ  (5, 206)  DIAM , DEPTH, TSOD 
READ  (5,20  6)  Q ,TSD,PSD,TDD 
IE(ITYPE)  1,1,2 


\  c  I 
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1  WRITE  (6,218) 

GC  TO  3 

2  WRITE (6,219) 

WRITE  (6,20  8) 

3  WRITE  (6 ,2 1 6)  BIAM, DEPTH 

W RITE  (6,217)  Q,TSOD,TSD, TDD , PSD 
CALL  AMEIN 

WRITE  (6,2  27)  SG , PTEMC , PPREC 
WRITE  (6,225) 

WRITE  (6,226) 

C  CALCULATING  MOLAL  FLOW  RATE 

TR=(0.6D  C 2+0 •  45967D  03)/ZTEMC 
P  R  =  0  . 1465D  0 2/ZPREC 

CALL  BWR  (TR, PR, ITER, ZST , DR , C PD M , DZDTR , DZDPR) 
VMASS=0. 55744134D-02*Q/  (ZST*DIAM**2) 


C  READ  1  CARD  (  RANGE  OF  R**  ) 

C  XMIN=1000.C 

C  XOUT= 1010.0 

C  XMAX=7000 . C 

C  IF  CONVERGENCE  IS  NOT  ACHIEVED, THE  RANGE 

C  MAY  BE  PUT  BETWEEN  200.0  &  1000.0  OR 

C  70C0.0  &  2C0 000.0 

**  REMARK  **  FOR  SAVING  COMPUTER  TIME  IN  CASE  R** 

IS  OUTSIDE  THE  SPECIFIED  RANGE, 

ALWAYS  PUT  XOUT=XMIN+1C. 0 

RE A D  (5 , 20 1 )  XMIN ,XOUT,XMAX 
TEMESP=C. ID  00 

C  MINIMIZING  FUNCTION  F 

CALL  SEARCH  ( F , TEMESP , X , XMI N , XM A X , XOUT) 

201  FORMAT (3F10. 4) 

2C2  FORMAT  (10X,3F16. 7) 

204  FORMAT (4D15.7) 

205  FCRMAT  (10X,4D20.7) 

206  FORMAT  (5F15. 8) 

2C7  FORMAT  (18X,F20. 5) 

208  FORMAT  (' 1 ') 

209  FCRMAT  (415, 3F10. 4) 

211  FORMAT (10X,5F15. 7) 

213  FORMAT (///////, 29X , • A ■ , 1 9X, • B ' , 1 9X, ■ C • , 1 9X, • D ' , /) 

214  FCRMAT (//,47X,'CPO (I) =A+B*T+C*T**2+D*T**3 • ) 

215  FORMAT  (' 1 ',////////, 27X, • MOLE  FRACTION',/) 

216  FORMAT  (//////, 10X, 'TUBING  I.  D. ' , 27X , '  =  ' , FI  0 . 3 , 

A  '  INCHES',/, 

B1 OX , ' THE  LENGTH  OF  FLOW  STRING  =', 


(N  CN  CM 


C  FI 0.3,'  FEET') 

217  FORMAT  (1  OX,*  FLOW  RATE ' , 30X , '  =  ' ,F 10. 3 , 

A  '  MSCF/DAY ' ,/, 

B10X, 'THE  AVE.  TEMP.  OF  THE  GROUND  SURFACE  =', 

C  F10.3, '  F ' , /, 

D10X, 'THE  TEMP.  OF  FLOWING  STREAM  AT  SURFACE  =', 

E  F10.3, *  F ' , / , 

F 1 0X  ,  'THE  TEMP.  OF  THE  PRODUCING  FORMATION  =', 

G  F 1  0 . 3  ,  '  F  '  ,  /  , 

H 1 0X , ' THE  WELL  HEAD  TUBING  PRESSURE  =*, 

I  F10 . 3 , *  PSIA',/) 

18  FORMAT(10X, 'INJECTION  PROBLEM') 

19  FORMAT (1CX, ' PRODUCTION  PROBLEM') 

25  FORMAT  (10X, 'THE  CALCULATED  BOTTOM-HOLE  PRESSURE' 

1  *  IS  THAT  WHEN  THE  ABSOLUTE  DIFFERENCE  BETWEEN' 

2  /,10X, 'EARTH  AND  CALCULATED  FLUID  TEMPERATURES' 

3  '  AT  THE  BOTTOM  OF  WELL  IS  WITHIN  0.1  DEG.  F  .' 

226  FORMAT (/, 10X, 'FIRST  COLUMN  IS  DEPTH  IN  FEET  .',/ 

1  10X, 'SECOND  COL.  IS  EARTH  TEMP.  IN  DEG.  F  .',/, 

2  10X, 'THIRD  COL.  IS  CALC.  FLUID  TEMP.  IN  DEG.  F 

3  ,/,10X, ' FOURTH  COL.  IS  CALC.  PRESS.  IN  PSIA  .', 

4  /,10X, 'FIFTH  COL.  IS  R**  .',/) 

227  FORMAT  (10X, 'GAS  GRAVITY,  PSEUDO  REDUCED  TEMP.,', 
1  'AND  PSEUDO  REDUCED  PRESS.  ARE  : ' , /, 1 0 X , 3F1 0 . 4 ) 

STOP 

END 


»%e.ors  o 

: 

%\% •yAa\^D2K 

%*  =  30A9HUE  TA  MA39T2  0HIHOI3  iO  .9M3I  3HT'»X0ra 

*%e.on  3 

OI  :  •:  H 1003039  I  r.O  .'if..  :  IB  »%XOrr? 

At*?  ’A.Or*  D 

,»=  23032393  OKU,.'  GA'H  IJ3W  9flT',X0rH 

•A. or*  i 

(•81318099  WO”  I ;.  •  r )  TAM903  8T£ 

(  :  .  a  jh‘-  -  v  ,••  ”■  r )  r  •  o  i  ers 

%  »30O827£t:3  3,1'  I  - iOTTOQ  dSr  ?  (JZ  ;  HI'*  4X0  r  VAbHOl  iZS. 
,»K;[2WT20  9011333131(1  3TUI02HA  SHf  H;  -i r-'  T  , l  :  21  •  T 
,  •  3  ?  i  OTA 1  '.3  9  M  ST  a  1 0  IS  3  ST  A I UOI  i  0  Q  ;  H  T9  A  2  •  ,  7  0  T  .  ,  2 
(«  .  3  .OSQ  r.O  MIHTX 21  IIS  10  ’TOL  3HT  TA  •  Z 
A%  1  •  T333  HI  hT'3’13  21  HMUJOO  rl  '  V  %\)  TASSO  I  d££ 
,  \  »  *  •  3  .osa  K  .  a:  U  HT-A2  U  .  I  >,*  Qli  0032  » -  X  AT  f 
•.  ri  .osa  HI  .9 M3T  (II0I3  .01  AO  21  .100  1  :  i  T  *  „  ,  '  £ 

,•*  AI39  KI  .22393  .OIAD  21  .100  hTBUOS  •  ,  XOI\\%  £ 

(\,».  **9  3i  .loo  mi  ,xOr,\  # 

,  •  %  .9MST  asouaas  oausas  ,ytu  u  3Ai  . or  tahho  i  rss 

(^.Or3e,xorA%* :  ssa  .33399  aaonuu  00  .  q  ika*  r 


90T3 

CIKS 


160 


FUNCTION  F  (X) 

C****  THIS  FUNCTION  SUBROUTINE  IS  TO  BE 
C****  USED  WITh  MAIN  MODEL  BB  ONLY 

C  F  IS  THE  FUNCTION  TO  BE  MINIMIZED 

C  E  IS  ABSOLUTE  ROUGHNESS  OF  TUBING  IN  INCHES 

C  DIAM  IS  TUBING  DIAMETER  IN  INCHES 

C  DL  IS  MONITORED  INTEGRATION  STEP  SIZE  IN  FEET 

C  TGEOG  IS  GEOTHERMAL  GRADIENT  ,  DEG.  F/  FOOT 

C  PRE  IS  PRESSURE  IN  PSIA 

C  TEM  IS  TEMPERATURE  IN  DEG.  RANKIN 

C  DENS  IS  DENSITY  OF  NAT.  GAS  IN  GM/C.C. 

C  REN  IS  REYNOLDS  NUMBER 

C  OHTCO  IS  THE  OVERALL  HEAT-TRANSFER  COEFFICIENT 

IMPLICIT  REAL *8  (A-H,0-Z) 

DIMENSION  PRY  (5)  ,  TRY  (5)  ,PRK  (5)  ,TRK  (5) 

COMMON  /BL3/AVEMW/BL1 3/DIAM/BL1 4/SG/BL1 8/M 

COMMON  /BL12/PTEMC, ZTEMC , PPREC, Z PREC 

COMMON  /BL 1 5/Q ,TSD, PSD, TDD, DEPTH, TSOD,DL 

COMMON  /BL30/VMASS 

COMMON  /BL50/MMM 

D1=0. 1D-0  1 

D2=0. ID-01 

DL=C • 1 D  02 

E=G . 6D-03 

RST AR=X 

RR=E/DIAM 

TGEOG= (TDD-TSOD) /DEPTH 

C  THE  RUNGE-KUTTA-MERSON  INTEGRATION  ALGORITHM 

TEM=TSD+0. 45967D  03 
PRE=PSD 
T  S=T  SOD 
TRN=TEM/PTEMC 
PRN= ERE/P PREC 
TL=C.0D  CC 

W  RITE (6,211)  TL,TS,TSD,PRE,RSTAR 

1  DO  14  Jl  =  1 ,5 

IF  (Jl-1)  2,2,3 

2  PRY  (1)  =PR N 
TRY  (1)  =TRN 

TS=T  SOD+TGEOG*TL 

3  IF  ( Jl-2)  11,4,5 

4  PRY  (2)  =PRN  +  PRK  (1) 

TRY  (2)  =TRN+TRK  (1) 

TS=TSOD+TGEOG* (TL+DL/0.  3D  01) 

5  IF  ( Jl-3)  11,6,7 

6  PRY  (3)  =PRN+0.  5D  00*PRK  (1 ) +0 . 5D  0Q*PRK(2) 

TRY  (3)  =TRN  +  0.  5D  00*TRK ( 1 ) +0 . 5D  00*TRK(2) 
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TS=TSOD+TGEOG* (TL+DL/O. 3D  Cl) 

7  IF  (J1-4)  11,8,9 

8  PRY  (4)  =PEN+0. 375D  00*PRK  (1 ) +0 . 1 1  25D  01*PRK(3) 
TRY  (4) =TRN+0. 375D  00*TRK ( 1 ) +0 . 1 1 25D  01*TRK(3) 
TS=TSOD+TGEOG*  (TL+DL/0.  2D  01) 

9  IF  (J1-5)  11,10,14 

10  ABC1=PRN+C. 15D  0 1  *PRK  ( 1 ) 

ABC2=-0.45D  0 1*PRK  (3) +0. 6D  01*PRK(4) 

PRY  (5) = ABC1 +ABC2 

ABC1=TRN+C. 15D  01*TRK(1) 

AEC2=-G.45D  0 1*TRK  (3) +0. 6D  01*TRK(4) 

TRY  (5)  =ABC1+ABC2 
T  S  =  T SOD  +  TGEOG*  (TL+DL) 

PRJ AM=PRY  (5) 

TR J AM=TRY  (5) 

11  PR=PRY(J1) 

TR  =  TRY  (J1 ) 

TRS=  (TS  +  G.  45967D  03)/PTEMC 

PRE=PR*PPREC 

TEM=TR*PTEMC 

C  CALCULATING  IDEAL  HEAT  CAPACITY  OF  NAT.  GAS 

CALL  KASEM  (TR,TEM,SG,M,CPOM) 

TRC=TEM/ZTEMC 

PRC=PRE/ZPREC 

C  CALCULATING  COMPRESSIBILITY  FACTOR, ITS  PARTIAL 

C  DERIVATIVES  W.R.T.  PR  &  TR  , REDUCED 

C  DENSITY,  AND  HEAT  CAPACITY  DEPARTURE 

C  CF  NAT.  GAS 

CALL  BWR  (TRC , PRC , ITER, Z , DR ,CPDM ,DZDTR, DZDPR) 
IF  (ITER)  22,22,12 

12  DEN S=0 •  14925826D-Q2*AVEMW*PEE/(Z*TEM) 
CPM=CPDM+CPOM 

C  CALCULATING  VISCOSITY  OF  NAT.  GAS  IN  C.  P. 

CALL  LEE  (TEM,AVEMW,DENS, VISM) 

REN  =  0. 200 1485 5D  02*Q*SG/  (VISM*DIAM) 

C  CALCULATING  THE  HYDRAULIC  FRICTION  FACTOR 

CALL  MOODY  (RR , REN , FF , I) 

U 1 =RSTAR 

OHTCO=0 • ID  01/U1 

SOLVING  FOR  THE  PATE  OF  CHANGE  OF  REDUCED 
PRESSURE  AND  TEMPERATURE  WITH  LENGTH 
CALL  QUON  (TEM,PRE,TR, PR, TRS,Z, DZDPR, DZDTR, 
lFF,CPM,OHTCO,VMASS,DPRDL,DTRDL, NASER) 
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IF(NASER-I)  22,13,13 

13  PRK  (J1) =DL*DPRDL/0. 3D  01 
TRK (J1) =DL*DTRDL/0. 3D  01 

14  CONTINUE 

MONITORING  THE  INTEGRATION  STEP  SIZE 

ABC5=PRK  (1) -0. 45D  01*PRK(3) 

A BC6=C • 4 D  01 *PRK  (4)  -0 . 5D  00*PRK  (5) 

EPR  =  ABC  5  +  ABC6 

ABC5=TRK  (1) -0. 45D  01*TRK(3) 

ABC6=0.4D  01*TRK  (4) -0.5D  CC’*TRK  (5) 

ETR= ABC5+ ABC6 
ETEM=DABS  (ETR*PTEMC) 

EPRE  =  DABS  (EPR*PPREC) 

IF  (EPRE-D  1 )  15,15,16 

15  IF  (ETEM-D2)  17,17,16 

16  DL=DL/0.2D  01 
GO  TO  1 

17  P R=PRJ AM- EPR 
TR=TRJAM-ETR 
TI=TL+DL 

IF (EPRE-D1/Q. 16D  02)  18,20,20 

18  IF  (ETEM-D2/0. 16D  02)  19,20,20 

19  DL=0. 2D  0 1 *DL 
2C  PRN=PR 

TRN=TR 

PRE=PRN*PPREC 

TEM=TRN*PTEMC 

TED=TEM-0. 45967D  03 

WRITE  (6,211)  TL,TS,TED,PRE,RSTAR 

DEI=DEPTH-TL 

IF  (CABS  (DEL) -0. ID-04)  23,23,25 
25  I F  ( CEL-DL )  21,1,1 

21  DL=DEL 
GO  TO  1 

22  WRITE  (6,2  27) 

23  TED=TEM-C • 45967D  03 
TEMTL=TED 
F=TED-TEMTL 

211  FORMAT  (10X,5F10. 2) 

227  FORMAT (//, 10X  ,  ' ***  THE  COURSE  OF  ITERATIONS  IS 
INTERRUPTED  ****,//, 7X ***  THE  PROGRAM  IS' 
*,»  GOING  TO  ASSUME  ANOTHER  "RSTAR"*** ' ) 

RETURN 

END 
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SUBROUTINE  AMEIN 

C****  THIS  SUBROUTINE  CAN  BE  USED  WITH 
C****  MAIN  MODEL  BB  ONLY 


THIS  SUBROUTINE  CALCULATES  THE  MOLAL  AVERAGE 

PROPERTIES  AND  THE  WICHERT  &  AZIZ  CORRECTED 
PSEUDCCRITICAL  PROPERTIES  OF  NAT.  GAS 


AVEMW  IS 
SG  IS 
PPREC  IS 
PTEMC  IS 
ZTEMC  IS 
ZPREC  IS 


AVERAGE  MOLECULAR  WEIGHT 
GAS  GRAVITY 

PSEUDOCRITICCAL  PRESSURE  IN  PSIA 
PSEUDOCRITICAL  TEMPERATURE , DEG .  R 
CORRECTED  PSEUDOCRITICAL  TEMPERATURE 
CORRECTED  PSEUDOCRITICAL  PRESSURE 


IMPLICIT  REAL*8  (A-H,0-Z) 

COMMON  /BL3/A VEM W/BL5/YM W (15) /BL9/YMOLF (1 5) 
COMMON  /3L1 1/TEMC  (15)  ,PREC  (15) 

COMMON  /BL12/PTEMC, ZTEMC, PPREC, ZPREC 
COMMON  /BL14/SG/BL1 8/M/BL3 1/MM 
COMMON  /BL50/MMM 
Y 1 1 =YMCLF  (1) 

Y22=YMOLF  (2) 

Y33=YMOLF  (3) 

IF(MMM-I)  10,3,3 
1C  IF(M-I)  3,1,1 

1  AVEMW=0.0D  00 
PTE MC=C .CD  00 
PPREC=0. CD  00 
DO  2  1=1,15 

A VEMW  =  A VEMW+ YMOLF (I)  *YMW (I) 

PTEMC=PTEMC  +  YMOLF  (I) *TEMC (I) 

2  PPREC=PPR EC+YMOLF  (I) *PREC  (I) 

SG=AVEMW/0 . 2397D  02 

GO  TO  7 

3  AVEMW=SG*C.2897D  02 

IF  (MM- 1 )  11,7,7 

THE  BROWN  ET  AL  AND  CARR  ET  AL  CORRELATIONS 
FOR  THE  PSEUDOCRITICAL  PROPERTIES  OF 
SWEET  AND  SOUR  NAT.  GASES 


11  PTEMC=0. 1 71137D  03+0.313725D  03*SG 

TCCOR=-0.8D  Q2*Y11-0.25D  03*Y22+0.13D  03*Y33 

PTEMC=PTEMC+TCCOR 

IF  (SG-0. 85D  00)  4,4,5 

4  PPREC=0 • 6 951 D  C3-C.4D  02*36 
GO  TO  6 

5  PPREC=0.7C4396D  03-0. 51724D  C2*SG 

6  PCCCR=0 . 4  4D  03*Y11-0.17D  03* Y22  +  C. 6D  C3*Y33 
PPREC=PPREC+PCCOR 
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C  THE  WICHE5T  &  AZIZ  CORRECTED  PSEUDOCRITICALS 

C  DUE  TO  PRESENCE  OF  C02  AND  H2S 

7  Sl=YMOLF ( 1) +YMOLF (3) 

S  2=  Y  MOLF (3) 

E31 =0. 1  2D  03*  (S1**0. 9D  0Q-S1**G.16D  01) 

E32  =  0.15D  02*  (DSQRT  (S2) -S2**4) 

E3=  E31 +  E3  2 

S4  =  FTEMC+S2*  (0, ID  01-S2)*E3 

ZTE  MC=PTEMC-E3 

ZPREC=PPR  EC*ZTEMC/S4 

RETURN 

END 


